Study overview: This study was conducted at two location along the Po river in Northwestern Italy. The two sites were chosen to represent two distinct riparian condition. One site (Pian della Regina) was an alpine prairie with very sparse vegetation located above the treeline. The second site (Ostana) consisted of a dense riparian forest and much higher levels of shade and leaf inputs. Macroinvertebrate samples were hand collected and identified to species as well as classified to macroinvertebrate functional feeding groups.
After the surfaces of the macroinvertebrates were decontaminated, the internal bacterial communities were identified throught high throughtput sequencing of the 16S V4 region (515f, 806r) on an Illumina MiSeq platform. After filtering in QIIME2 (2018.11), bacterial reads were classified to taxanomic and functional orthologs (KEGG) using QIIME 2 and PiCRUSt 2.
A total of 26 invertebrates were sequenced (see metadata tab below for breakdown)
To see the corresponding code for each output click on the button labeled code above each box. All error bars are SEM and all multiple comparisons are adjusted with a FDR correction unless otherwise noted.
For code used to make figures see “ManuscriptFigures.r”. The .biom files and metadata are in the Github project repository.
Data overview
Sequencing
Based on the alpha rarefaction curves (Figure S1), samples were rarefied to 3,000 reads to limit the impact of differential library sizes on diversity metrics.
physeq## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2420 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 2420 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 2420 tips and 2400 internal nodes ]Metadata
On average, predators, scrapers, and shredders were larger at the forested station (Ostana) than at the alpine prairie site.
Trtdata <- ddply(metadata, c("Sampling_station", "FFG","Taxon_name"), summarise,
N = length(id),
meanWeight = mean(Mass),
sd = sd(Mass),
se = sd / sqrt(N)
)
Trtdata## Sampling_station FFG Taxon_name N meanWeight sd
## 1 Ostana Filterer Hydropsyche 4 41.4500 33.707615
## 2 Ostana Predator Perla 2 190.5000 174.513954
## 3 Ostana Scraper Epeorus 5 22.7600 3.857849
## 4 Ostana Shredder Tipula_(Acutitipula) 3 534.9333 425.432889
## 5 PDR Predator Dictyogenus 4 83.9500 16.270730
## 6 PDR Scraper Epeorus 4 11.6000 4.009988
## 7 PDR Shredder Tipula_(Acutitipula) 4 87.7000 42.817987
## se
## 1 16.853808
## 2 123.400000
## 3 1.725283
## 4 245.623793
## 5 8.135365
## 6 2.004994
## 7 21.408993ggplot(metadata,aes(x=Sampling_station,y=Mass,fill=Sampling_station))+geom_boxplot()+ #Write in labels from posthoc
scale_fill_manual(values=cbPalette)+xlab("")+ylab("Mass (mg)")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG,scales = "free_y")#+theme(axis.text.x = element_text(angle = 45, hjust = 1))+#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
head(metadata)## id FieldID Taxon_name FFG Sampling_station Mass
## 1 II5 T1OS Tipula_(Acutitipula) Shredder Ostana 989.1
## 2 II6 T2OS Tipula_(Acutitipula) Shredder Ostana 470.0
## 3 II8 T4OS Tipula_(Acutitipula) Shredder Ostana 145.7
## 4 II13 P1OS Perla Predator Ostana 67.1
## 5 II14 P2OS Perla Predator Ostana 313.9
## 6 II15 H2OS Hydropsyche Filterer Ostana 7.2
## shannon faith_pd
## 1 7.086586 20.95169
## 2 7.430371 22.58281
## 3 7.383975 21.47860
## 4 6.952585 17.94331
## 5 6.200620 26.41718
## 6 6.153952 17.75811metadataStatsSubset<- subset(metadata,FFG!="Predator")
metadataStatsSubset<- subset(metadataStatsSubset,FFG!="Filterer")
metadataStatsSubset## id FieldID Taxon_name FFG Sampling_station Mass
## 1 II5 T1OS Tipula_(Acutitipula) Shredder Ostana 989.1
## 2 II6 T2OS Tipula_(Acutitipula) Shredder Ostana 470.0
## 3 II8 T4OS Tipula_(Acutitipula) Shredder Ostana 145.7
## 10 E1OS E1OS Epeorus Scraper Ostana 20.2
## 11 E2OS E2OS Epeorus Scraper Ostana 27.5
## 12 E3OS E3OS Epeorus Scraper Ostana 19.0
## 13 E4OS E4OS Epeorus Scraper Ostana 20.8
## 14 E5OS E5OS Epeorus Scraper Ostana 26.3
## 15 II1 T1P Tipula_(Acutitipula) Shredder PDR 131.8
## 16 II2 T2P Tipula_(Acutitipula) Shredder PDR 53.9
## 17 II3 T3P Tipula_(Acutitipula) Shredder PDR 48.2
## 18 II4 T4P Tipula_(Acutitipula) Shredder PDR 116.9
## 23 E1P E1P Epeorus Scraper PDR 11.0
## 24 E3P E3P Epeorus Scraper PDR 17.4
## 25 E4P E4P Epeorus Scraper PDR 9.6
## 26 E5P E5P Epeorus Scraper PDR 8.4
## shannon faith_pd
## 1 7.0865864 20.951686
## 2 7.4303709 22.582812
## 3 7.3839746 21.478599
## 10 4.9093727 11.372657
## 11 1.9619653 9.905822
## 12 3.0828744 9.146140
## 13 4.8326128 7.641727
## 14 3.0611306 7.197870
## 15 7.2937707 21.728816
## 16 6.9992923 20.798101
## 17 6.0582633 11.170252
## 18 7.1744810 21.997409
## 23 1.6024907 4.266886
## 24 4.5283738 5.468116
## 25 0.7750815 4.605169
## 26 0.7526593 4.586619print("Mass (mg)~ Sampling_station*FFG, method=ANOVA")## [1] "Mass (mg)~ Sampling_station*FFG, method=ANOVA"av=aov( Mass~ Sampling_station*FFG, data=metadataStatsSubset)
summary(av)## Df Sum Sq Mean Sq F value Pr(>F)
## Sampling_station 1 109131 109131 3.563 0.08351 .
## FFG 1 319410 319410 10.427 0.00723 **
## Sampling_station:FFG 1 184026 184026 6.007 0.03054 *
## Residuals 12 367594 30633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("(Since the interaction between sampling station and FFG was signinficant below are t-tests split out by FFG/n Mass (mg) ~Sampling-station (grouped by FFG, FDR correction)")## [1] "(Since the interaction between sampling station and FFG was signinficant below are t-tests split out by FFG/n Mass (mg) ~Sampling-station (grouped by FFG, FDR correction)"Scrapers<-subset(metadataStatsSubset, FFG=="Scraper")
t.test( Mass~ Sampling_station, data=Scrapers)##
## Welch Two Sample t-test
##
## data: Mass by Sampling_station
## t = 4.2191, df = 6.4396, p-value = 0.004753
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.79325 17.52675
## sample estimates:
## mean in group Ostana mean in group PDR
## 22.76 11.60Shredders<-subset(metadataStatsSubset, FFG=="Shredder")
t.test( Mass~ Sampling_station, data=Shredders)##
## Welch Two Sample t-test
##
## data: Mass by Sampling_station
## t = 1.8139, df = 2.0304, p-value = 0.2095
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -598.5184 1492.9851
## sample estimates:
## mean in group Ostana mean in group PDR
## 534.9333 87.7000#compare_means(Mass_.mg. ~ Sampling_station,group.by = "FFG", data = metadata, p.adjust.method = "fdr",method="t.test")Alpha Diversity
All three diversity metrics showed a similar pattern with scrapers having lower diversity than the other three groups.
Observed Species
plot_richness(physeq, x="FFG",color="FFG", measures=c("Observed"))+geom_boxplot(aes(x=FFG, y=value, color=FFG), alpha=0.05)+ylab("Observed Species")+geom_point(size=3)+guides(fill=FALSE)
Shannon Diversity
av=kruskal.test( shannon~ FFG, data=metadataStatsSubset) #fig.width = 8, fig.height = 8 for line above to change size
av2=kruskal.test( shannon~ Sampling_station, data=metadataStatsSubset)
av##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 11.118, df = 1, p-value = 0.0008551av2##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 0.54044, df = 1, p-value = 0.4622posthoc<-compare_means(shannon ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(metadataStatsSubset,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
stats2<-summarySE(metadata,measurevar="shannon",groupvars=c("FFG","Sampling_station",na.rm=TRUE))
posthoc2<-compare_means(shannon ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr",group.by = "Sampling_station")
ggplot(metadata,aes(x=FFG,y=shannon,fill=FFG))+geom_boxplot() +#Write in labels from posthoc+geom_col()+geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black")+
scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~Sampling_station)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
ggplot(metadata,aes(x=Sampling_station,y=shannon,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
av##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 11.118, df = 1, p-value = 0.0008551posthoc## # A tibble: 1 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 shannon Scraper Shredder 0.00103 0.001 0.001 ** WilcoxonLetters## Scraper Shredder
## "a" "b"stats## FFG na.rm N shannon sd se ci
## 1 Scraper NA 9 2.834062 1.6633640 0.5544547 1.2785748
## 2 Shredder TRUE 7 7.060963 0.4686418 0.1771300 0.4334214posthoc2## # A tibble: 2 x 9
## Sampling_station .y. group1 group2 p p.adj p.format p.signif
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 Ostana shan~ Scrap~ Shred~ 0.0369 0.037 0.037 *
## 2 PDR shan~ Scrap~ Shred~ 0.0304 0.037 0.030 *
## # ... with 1 more variable: method <chr>stats2## FFG Sampling_station na.rm N shannon sd se ci
## 1 Filterer Ostana NA 4 6.502903 0.3062089 0.1531045 0.4872467
## 2 Predator Ostana NA 2 6.576603 0.5317195 0.3759825 4.7773101
## 3 Predator PDR NA 4 6.164048 1.3841309 0.6920654 2.2024611
## 4 Scraper Ostana NA 5 3.569591 1.2718220 0.5687761 1.5791756
## 5 Scraper PDR NA 4 1.914651 1.7867880 0.8933940 2.8431785
## 6 Shredder Ostana TRUE 3 7.300311 0.1865387 0.1076982 0.4633878
## 7 Shredder PDR NA 4 6.881452 0.5619605 0.2809802 0.8942045av2##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 0.54044, df = 1, p-value = 0.4622AlphaDivSplitByFFG
print("Scrapers")## [1] "Scrapers"Scrapers<-subset(metadataStatsSubset, FFG=="Scraper")
shapiro.test(Scrapers$shannon)##
## Shapiro-Wilk normality test
##
## data: Scrapers$shannon
## W = 0.8943, p-value = 0.2208shapiro.test(Scrapers$faith_pd)##
## Shapiro-Wilk normality test
##
## data: Scrapers$faith_pd
## W = 0.91387, p-value = 0.3439t.test( shannon~ Sampling_station, data=Scrapers)##
## Welch Two Sample t-test
##
## data: shannon by Sampling_station
## t = 1.5626, df = 5.2748, p-value = 0.1759
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.025406 4.335285
## sample estimates:
## mean in group Ostana mean in group PDR
## 3.569591 1.914651t.test( faith_pd~ Sampling_station, data=Scrapers)##
## Welch Two Sample t-test
##
## data: faith_pd by Sampling_station
## t = 5.3855, df = 4.8851, p-value = 0.003194
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.243935 6.398356
## sample estimates:
## mean in group Ostana mean in group PDR
## 9.052843 4.731698kruskal.test( shannon~ Sampling_station, data=Scrapers)##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 2.94, df = 1, p-value = 0.08641kruskal.test( faith_pd~ Sampling_station, data=Scrapers)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 6, df = 1, p-value = 0.01431print("Shredders")## [1] "Shredders"Shredders<-subset(metadataStatsSubset, FFG=="Shredder")
t.test(shannon~Sampling_station,data=Shredders)##
## Welch Two Sample t-test
##
## data: shannon by Sampling_station
## t = 1.392, df = 3.8225, p-value = 0.2394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4321367 1.2698544
## sample estimates:
## mean in group Ostana mean in group PDR
## 7.300311 6.881452t.test( faith_pd~ Sampling_station, data=Shredders)##
## Welch Two Sample t-test
##
## data: faith_pd by Sampling_station
## t = 1.0402, df = 3.2033, p-value = 0.3703
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.364469 10.859245
## sample estimates:
## mean in group Ostana mean in group PDR
## 21.67103 18.92364kruskal.test(shannon~Sampling_station,data=Shredders)##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 2, df = 1, p-value = 0.1573kruskal.test( faith_pd~ Sampling_station, data=Shredders)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.5, df = 1, p-value = 0.4795AlphaDivSplitByStation
print("Ostana")## [1] "Ostana"Ostana<-subset(metadata, Sampling_station=="Ostana")
OstanaStatAlpha<-subset(Ostana,FFG!="Predator")
shapiro.test(OstanaStatAlpha$shannon)##
## Shapiro-Wilk normality test
##
## data: OstanaStatAlpha$shannon
## W = 0.87409, p-value = 0.07366shapiro.test(OstanaStatAlpha$faith_pd)##
## Shapiro-Wilk normality test
##
## data: OstanaStatAlpha$faith_pd
## W = 0.8379, p-value = 0.02611kruskal.test( shannon~ FFG, data=OstanaStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 9.6923, df = 2, p-value = 0.007859compare_means(shannon ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 shannon Filterer Scraper 0.0200 0.052 0.020 * Wilcoxon
## 2 shannon Filterer Shredder 0.0518 0.052 0.052 ns Wilcoxon
## 3 shannon Scraper Shredder 0.0369 0.052 0.037 * Wilcoxonkruskal.test( faith_pd~ FFG, data=OstanaStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 8.2564, df = 2, p-value = 0.01611compare_means(faith_pd ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 faith_pd Filterer Scraper 0.0200 0.055 0.020 * Wilcoxon
## 2 faith_pd Filterer Shredder 0.596 0.6 0.596 ns Wilcoxon
## 3 faith_pd Scraper Shredder 0.0369 0.055 0.037 * Wilcoxonprint("PDR")## [1] "PDR"PDRStatAlpha<-subset(metadata, Sampling_station!="Ostana")
kruskal.test( shannon~ FFG, data=PDRStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 7.7308, df = 2, p-value = 0.02095compare_means(shannon ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 shannon Predator Scraper 0.0304 0.046 0.03 * Wilcoxon
## 2 shannon Predator Shredder 0.470 0.47 0.47 ns Wilcoxon
## 3 shannon Scraper Shredder 0.0304 0.046 0.03 * Wilcoxonkruskal.test( faith_pd~ FFG, data=PDRStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 7.5385, df = 2, p-value = 0.02307compare_means(faith_pd ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 faith_pd Predator Scraper 0.0304 0.046 0.03 * Wilcoxon
## 2 faith_pd Predator Shredder 0.665 0.67 0.67 ns Wilcoxon
## 3 faith_pd Scraper Shredder 0.0304 0.046 0.03 * Wilcoxonposthoc<-compare_means(shannon ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(OstanaStatAlpha,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity (Ostana)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
posthoc<-compare_means(shannon ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(PDRStatAlpha,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity (PDR)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
posthoc<-compare_means(faith_pd ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(OstanaStatAlpha,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=faith_pd,fill=FFG))+geom_text(aes(x=FFG, y=faith_pd+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=faith_pd-se,ymax=faith_pd+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD (Ostana)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
posthoc<-compare_means(faith_pd ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(PDRStatAlpha,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=faith_pd,fill=FFG))+geom_text(aes(x=FFG, y=faith_pd+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=faith_pd-se,ymax=faith_pd+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD (PDR)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
Faith’s PD
av=kruskal.test( faith_pd~ FFG, data=metadataStatsSubset)
av2=kruskal.test( faith_pd~ Sampling_station, data=metadataStatsSubset)
av##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 10.423, df = 1, p-value = 0.001245av2##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.89338, df = 1, p-value = 0.3446posthoc<-compare_means(faith_pd ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(metadata,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(metadata,aes(x=FFG,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD")+labs(fill="FFG")#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))geom_text(aes(x=FFG, y=faith_pd+se+2,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ 
stats2<-summarySE(metadata,measurevar="faith_pd",groupvars=c("FFG","Sampling_station",na.rm=TRUE))
posthoc2<-compare_means(faith_pd ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr",group.by = "Sampling_station")
ggplot(metadata,aes(x=FFG,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith PD")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~Sampling_station,scales = "free_x")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
ggplot(metadata,aes(x=Sampling_station,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
av##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 10.423, df = 1, p-value = 0.001245posthoc## # A tibble: 1 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 faith_pd Scraper Shredder 0.00150 0.0015 0.0015 ** WilcoxonLetters## Scraper Shredder
## "a" "b"stats## FFG na.rm N faith_pd sd se ci
## 1 Filterer NA 4 21.911376 2.894769 1.4473845 4.606223
## 2 Predator NA 6 21.557631 5.684157 2.3205474 5.965157
## 3 Scraper NA 9 7.132334 2.594210 0.8647367 1.994086
## 4 Shredder TRUE 7 20.101096 3.984788 1.5061083 3.685314posthoc2## # A tibble: 2 x 9
## Sampling_station .y. group1 group2 p p.adj p.format p.signif
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 Ostana fait~ Scrap~ Shred~ 0.0369 0.037 0.037 *
## 2 PDR fait~ Scrap~ Shred~ 0.0304 0.037 0.030 *
## # ... with 1 more variable: method <chr>av2##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.89338, df = 1, p-value = 0.3446Taxa Plots
GPr = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPr = filter_taxa(GPr, function(x) mean(x) > 1e-5, TRUE)
PhylumAll=tax_glom(GPr, "Phylum")
PhylumLevel = filter_taxa(PhylumAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower tha 0.1%
FamilyAll=tax_glom(GPr,"Family")
FamilyLevel = filter_taxa(FamilyAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%
GenusAll=tax_glom(GPr,"Genus")
GenusLevel = filter_taxa(GenusAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%Phylum Level
The top 7 phyla across all samples are given in the table below
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
df2 <- psmelt(PhylumAll)
df2$Abundance<-df2$Abundance*100
PhylumAll<-ddply(df2, c("Phylum"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
PhylumSorted<-PhylumAll[order(-PhylumAll$mean),]
PhylumSorted[1:7,]## Phylum N mean sd se
## 20 Proteobacteria 26 51.896154 21.817321 4.2787287
## 5 Bacteroidetes 26 17.466667 14.180099 2.7809462
## 13 Firmicutes 26 13.275641 17.088823 3.3513939
## 23 Tenericutes 26 10.424359 19.936032 3.9097776
## 26 Verrucomicrobia 26 1.524359 1.846861 0.3621992
## 19 Planctomycetes 26 1.494872 1.825686 0.3580464
## 3 Actinobacteria 26 1.121795 1.100947 0.2159134PhylumAll2<-ddply(df2, c("Phylum","FFG"), summarise, #For Rel abu tabs
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
#PhylumAll2 Run this to see all taxa, commented out since its a large listcdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)
cdataplot#+ scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(
Trtdata2 <- ddply(df, c("Phylum", "FFG","Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
dfStats<-subset(df, FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station,scales = "free_x")
cdataplot2
########
compare_means(Abundance ~ FFG, data = df, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.0147 0.021 0.01466 * Kruskal-Wall~
## 2 Tenericutes Abundance 0.0643 0.075 0.06432 ns Kruskal-Wall~
## 3 Bacteroidetes Abundance 0.000400 0.0014 0.00040 *** Kruskal-Wall~
## 4 Firmicutes Abundance 0.000209 0.0014 0.00021 *** Kruskal-Wall~
## 5 Planctomycetes Abundance 0.0127 0.021 0.01274 * Kruskal-Wall~
## 6 Verrucomicrobia Abundance 0.00184 0.0043 0.00184 ** Kruskal-Wall~
## 7 Actinobacteria Abundance 0.118 0.12 0.11762 ns Kruskal-Wall~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Phylum", p.adjust.method = "fdr")
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Phylum))
for (i in levels(Means$Phylum)){
Tax<-i
TaxAbundance<-subset(Means,Phylum==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+10,label=vec)) +scale_fill_manual(values=cbPalette)
cdataplot
head(Trtdata)## Phylum FFG N mean sd se
## 1 Actinobacteria Filterer 4 0.3833333 0.4811252 0.2405626
## 2 Actinobacteria Predator 6 1.3111111 0.3862450 0.1576839
## 3 Actinobacteria Scraper 9 1.2259259 1.2891403 0.4297134
## 4 Actinobacteria Shredder 7 1.2476190 1.4698018 0.5555329
## 5 Bacteroidetes Filterer 4 22.6166667 3.7669616 1.8834808
## 6 Bacteroidetes Predator 6 32.7000000 10.7146006 4.3742174dfStatsPlot<-subset(Trtdata,Phylum!="Actinobacteria"&Phylum!= "Verrucomicrobia")
#vec<-c("a","b","a","b","a","b","a","b","a","b")
dfStatsPlot # Only show significant taxa ## Phylum FFG N mean sd se
## 5 Bacteroidetes Filterer 4 22.6166667 3.7669616 1.88348082
## 6 Bacteroidetes Predator 6 32.7000000 10.7146006 4.37421739
## 7 Bacteroidetes Scraper 9 0.8444444 0.6502136 0.21673788
## 8 Bacteroidetes Shredder 7 22.8380952 6.0532033 2.28789578
## 9 Firmicutes Filterer 4 20.0250000 13.0079285 6.50396425
## 10 Firmicutes Predator 6 0.2555556 0.4385160 0.17902341
## 11 Firmicutes Scraper 9 0.7888889 1.2328828 0.41096093
## 12 Firmicutes Shredder 7 36.6333333 7.7272486 2.92062543
## 13 Planctomycetes Filterer 4 1.3083333 0.4085884 0.20429418
## 14 Planctomycetes Predator 6 1.6000000 0.7636171 0.31174539
## 15 Planctomycetes Scraper 9 2.4444444 2.7655118 0.92183727
## 16 Planctomycetes Shredder 7 0.2904762 0.2052228 0.07756693
## 17 Proteobacteria Filterer 4 45.1250000 15.0889554 7.54447768
## 18 Proteobacteria Predator 6 51.1833333 9.3604072 3.82137022
## 19 Proteobacteria Scraper 9 68.0037037 27.1006890 9.03356301
## 20 Proteobacteria Shredder 7 35.6666667 9.4683488 3.57869948
## 21 Tenericutes Filterer 4 3.0750000 3.1347958 1.56739788
## 22 Tenericutes Predator 6 5.0666667 10.4961580 4.28503857
## 23 Tenericutes Scraper 9 24.9888889 28.2112665 9.40375549
## 24 Tenericutes Shredder 7 0.4904762 0.3207135 0.12121831cdataplot=ggplot(dfStatsPlot, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+#geom_text(aes(x=FFG, y=mean+se+10,label=vec))+
scale_fill_manual(values=cbPalette)+theme(legend.justification=c(0.05,0.95), legend.position=c(0.7,0.4))
cdataplot
#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Phylum))
#SigList<-length(unique(Trtdata$Phylum))
#SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
Means## # A tibble: 7 x 9
## Phylum .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobac~ Abunda~ Scraper Shred~ 0.0149 0.021 0.015 * Wilco~
## 2 Tenericut~ Abunda~ Scraper Shred~ 0.0148 0.021 0.015 * Wilco~
## 3 Firmicutes Abunda~ Scraper Shred~ 0.00103 0.0036 0.001 ** Wilco~
## 4 Bacteroid~ Abunda~ Scraper Shred~ 0.00103 0.0036 0.001 ** Wilco~
## 5 Planctomy~ Abunda~ Scraper Shred~ 0.0148 0.021 0.015 * Wilco~
## 6 Actinobac~ Abunda~ Scraper Shred~ 0.751 0.75 0.751 ns Wilco~
## 7 Verrucomi~ Abunda~ Scraper Shred~ 0.289 0.34 0.289 ns Wilco~MeansStation## # A tibble: 7 x 9
## Phylum .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacte~ Abunda~ Ostana PDR 0.958 1 0.958 ns Wilcox~
## 2 Tenericutes Abunda~ Ostana PDR 0.528 0.89 0.528 ns Wilcox~
## 3 Firmicutes Abunda~ Ostana PDR 1 1 1.000 ns Wilcox~
## 4 Bacteroidet~ Abunda~ Ostana PDR 0.637 0.89 0.637 ns Wilcox~
## 5 Planctomyce~ Abunda~ Ostana PDR 0.0180 0.13 0.018 * Wilcox~
## 6 Actinobacte~ Abunda~ Ostana PDR 0.495 0.89 0.495 ns Wilcox~
## 7 Verrucomicr~ Abunda~ Ostana PDR 0.0738 0.26 0.074 ns Wilcox~Means<-compare_means(Abundance ~ FFG, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
Trtdata <- ddply(dfStats, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Phylum))
for (i in levels(Means$Phylum)){
Tax<-i
TaxAbundance<-subset(Means,Phylum==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+5,label=vec))+scale_fill_manual(values=cbPalette)
cdataplot
#######Split by Station
Since there were not enough samples collect to compare between all groups at both stations, only the groups present with greater than 3 samples were compared at each station.
Ostana (Forest)
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")
dfStatsOstana<-subset(dfOstana,FFG!="Predator")
compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.291 0.290 0.2906 ns Kruskal-Wallis
## 2 Tenericutes Abundance 0.0549 0.077 0.0549 ns Kruskal-Wallis
## 3 Firmicutes Abundance 0.0121 0.042 0.0121 * Kruskal-Wallis
## 4 Bacteroidetes Abundance 0.00990 0.042 0.0099 ** Kruskal-Wallis
## 5 Planctomycetes Abundance 0.0430 0.077 0.0430 * Kruskal-Wallis
## 6 Actinobacteria Abundance 0.109 0.13 0.1089 ns Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.0484 0.077 0.0484 * Kruskal-WallisMeans<-compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## [1] Phylum group1 group2 p.adj
## <0 rows> (or 0-length row.names)Pian della Regina (Alpine Prairie)
print("Pian della Regina")## [1] "Pian della Regina"Means<-compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Phylum group1 group2 p.adj
## 2 Proteobacteria Predator Shredder 0.071
## 3 Proteobacteria Scraper Shredder 0.071
## 7 Bacteroidetes Predator Scraper 0.071
## 9 Bacteroidetes Scraper Shredder 0.071
## 11 Firmicutes Predator Shredder 0.071
## 12 Firmicutes Scraper Shredder 0.071
## 13 Verrucomicrobia Predator Scraper 0.071
## 14 Verrucomicrobia Predator Shredder 0.071
## 20 Planctomycetes Predator Shredder 0.071#########
############
ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Ostana(>1%)")
ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance PDR(>1%)")
print("Family level")## [1] "Family level"df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")
dfStatsOstana<-subset(dfOstana,FFG!="Predator")
print("Ostana")## [1] "Ostana"compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aeromonadaceae Abundance 0.0817 0.091 0.0817 ns Kruskal-Wal~
## 2 Comamonadaceae Abundance 0.00786 0.021 0.0079 ** Kruskal-Wal~
## 3 Rhodobacteraceae Abundance 0.886 0.89 0.8865 ns Kruskal-Wal~
## 4 Ruminococcaceae Abundance 0.00712 0.021 0.0071 ** Kruskal-Wal~
## 5 Desulfovibrionac~ Abundance 0.0105 0.021 0.0105 * Kruskal-Wal~
## 6 Lachnospiraceae Abundance 0.0118 0.021 0.0118 * Kruskal-Wal~
## 7 Rikenellaceae Abundance 0.0330 0.041 0.0330 * Kruskal-Wal~
## 8 Methylophilaceae Abundance 0.0139 0.021 0.0139 * Kruskal-Wal~
## 9 Enterobacteriace~ Abundance 0.0144 0.021 0.0144 * Kruskal-Wal~
## 10 Cytophagaceae Abundance 0.0139 0.021 0.0139 * Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Family group1 group2 p.adj
## 3 Comamonadaceae Filterer Scraper 0.058
## 4 Comamonadaceae Filterer Shredder 0.088
## 5 Comamonadaceae Scraper Shredder 0.076
## 9 Ruminococcaceae Filterer Scraper 0.054
## 11 Ruminococcaceae Scraper Shredder 0.054
## 12 Desulfovibrionaceae Filterer Scraper 0.054
## 14 Desulfovibrionaceae Scraper Shredder 0.054
## 15 Lachnospiraceae Filterer Scraper 0.054
## 17 Lachnospiraceae Scraper Shredder 0.054
## 18 Rikenellaceae Filterer Scraper 0.079
## 20 Rikenellaceae Scraper Shredder 0.054
## 21 Methylophilaceae Filterer Scraper 0.054
## 23 Methylophilaceae Scraper Shredder 0.063
## 25 Enterobacteriaceae Filterer Shredder 0.079
## 26 Enterobacteriaceae Scraper Shredder 0.076
## 27 Cytophagaceae Filterer Scraper 0.054
## 29 Cytophagaceae Scraper Shredder 0.063print("Pian della Regina")## [1] "Pian della Regina"compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteriace~ Abundance 0.0366 0.041 0.0366 * Kruskal-Wal~
## 2 Cytophagaceae Abundance 0.00609 0.013 0.0061 ** Kruskal-Wal~
## 3 Desulfovibrionac~ Abundance 0.00537 0.013 0.0054 ** Kruskal-Wal~
## 4 Lachnospiraceae Abundance 0.00921 0.013 0.0092 ** Kruskal-Wal~
## 5 Rikenellaceae Abundance 0.00537 0.013 0.0054 ** Kruskal-Wal~
## 6 Rhodobacteraceae Abundance 0.0249 0.031 0.0249 * Kruskal-Wal~
## 7 Comamonadaceae Abundance 0.00715 0.013 0.0072 ** Kruskal-Wal~
## 8 Ruminococcaceae Abundance 0.00921 0.013 0.0092 ** Kruskal-Wal~
## 9 Aeromonadaceae Abundance 0.368 0.37 0.3679 ns Kruskal-Wal~
## 10 Methylophilaceae Abundance 0.00609 0.013 0.0061 ** Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Family group1 group2 p.adj
## 3 Enterobacteriaceae Scraper Shredder 0.041
## 4 Cytophagaceae Predator Scraper 0.041
## 5 Cytophagaceae Predator Shredder 0.041
## 6 Cytophagaceae Scraper Shredder 0.041
## 7 Desulfovibrionaceae Predator Shredder 0.041
## 8 Desulfovibrionaceae Scraper Shredder 0.041
## 10 Lachnospiraceae Predator Shredder 0.041
## 11 Lachnospiraceae Scraper Shredder 0.041
## 12 Rikenellaceae Predator Shredder 0.041
## 13 Rikenellaceae Scraper Shredder 0.041
## 14 Rhodobacteraceae Predator Scraper 0.041
## 15 Rhodobacteraceae Predator Shredder 0.041
## 17 Comamonadaceae Predator Scraper 0.041
## 18 Comamonadaceae Predator Shredder 0.041
## 19 Comamonadaceae Scraper Shredder 0.041
## 21 Ruminococcaceae Predator Shredder 0.041
## 22 Ruminococcaceae Scraper Shredder 0.041
## 25 Methylophilaceae Predator Scraper 0.041
## 26 Methylophilaceae Predator Shredder 0.041
## 27 Methylophilaceae Scraper Shredder 0.041ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Ostana(>2%)")
ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance PDR(>2%)")
print("Genus Level")## [1] "Genus Level"df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")
dfStatsOstana<-subset(dfOstana,FFG!="Predator")
print("Ostana")## [1] "Ostana"compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 3 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacter Abundance 0.886 0.89 0.886 ns Kruskal-Wallis
## 2 Clostridium Abundance 0.0118 0.021 0.012 * Kruskal-Wallis
## 3 Methylotenera Abundance 0.0139 0.021 0.014 * Kruskal-Walliscompare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Genus", p.adjust.method = "fdr")## # A tibble: 9 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobact~ Abunda~ Filter~ Scraper 0.903 0.9 0.903 ns Wilcox~
## 2 Rhodobact~ Abunda~ Filter~ Shredd~ 0.860 0.9 0.860 ns Wilcox~
## 3 Rhodobact~ Abunda~ Scraper Shredd~ 0.766 0.9 0.766 ns Wilcox~
## 4 Clostridi~ Abunda~ Filter~ Scraper 0.0108 0.05 0.011 * Wilcox~
## 5 Clostridi~ Abunda~ Filter~ Shredd~ 0.596 0.89 0.596 ns Wilcox~
## 6 Clostridi~ Abunda~ Scraper Shredd~ 0.0168 0.05 0.017 * Wilcox~
## 7 Methylote~ Abunda~ Filter~ Scraper 0.0151 0.05 0.015 * Wilcox~
## 8 Methylote~ Abunda~ Filter~ Shredd~ 0.596 0.89 0.596 ns Wilcox~
## 9 Methylote~ Abunda~ Scraper Shredd~ 0.0262 0.059 0.026 * Wilcox~print("Pian della Regina")## [1] "Pian della Regina"compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 3 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Clostridium Abundance 0.00537 0.0091 0.0054 ** Kruskal-Wallis
## 2 Rhodobacter Abundance 0.0246 0.025 0.0246 * Kruskal-Wallis
## 3 Methylotenera Abundance 0.00609 0.0091 0.0061 ** Kruskal-Walliscompare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Genus", p.adjust.method = "fdr")## # A tibble: 8 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Clostridi~ Abunda~ Predat~ Shredd~ 0.0211 0.035 0.021 * Wilcox~
## 2 Clostridi~ Abunda~ Scraper Shredd~ 0.0211 0.035 0.021 * Wilcox~
## 3 Rhodobact~ Abunda~ Predat~ Scraper 0.0304 0.035 0.030 * Wilcox~
## 4 Rhodobact~ Abunda~ Predat~ Shredd~ 0.0294 0.035 0.029 * Wilcox~
## 5 Rhodobact~ Abunda~ Scraper Shredd~ 1 1 1.000 ns Wilcox~
## 6 Methylote~ Abunda~ Predat~ Scraper 0.0211 0.035 0.021 * Wilcox~
## 7 Methylote~ Abunda~ Predat~ Shredd~ 0.0304 0.035 0.030 * Wilcox~
## 8 Methylote~ Abunda~ Scraper Shredd~ 0.0211 0.035 0.021 * Wilcox~ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Genus)+ylab("Relative Abundance Ostana(>2%)")
ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Genus)+ylab("Relative Abundance PDR(>2%)")
Split by FFG
This set of tests is comparing bacterial between sites for a particular group (e.g. Shredders at Ostana va PDR)
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
dfShredders<-subset(df,FFG =="Shredder")
dfScrapers<-subset(df,FFG == "Scraper")
print("Shredders")## [1] "Shredders"compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.480 0.56 0.480 ns Kruskal-Wallis
## 2 Firmicutes Abundance 0.480 0.56 0.480 ns Kruskal-Wallis
## 3 Bacteroidetes Abundance 0.0339 0.18 0.034 * Kruskal-Wallis
## 4 Actinobacteria Abundance 0.0771 0.18 0.077 ns Kruskal-Wallis
## 5 Tenericutes Abundance 0.372 0.56 0.372 ns Kruskal-Wallis
## 6 Planctomycetes Abundance 0.0771 0.18 0.077 ns Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.589 0.59 0.589 ns Kruskal-WallisMeans<-compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Phylum group1 group2 p
## 3 Bacteroidetes Ostana PDR 0.05182993print("Scrapers")## [1] "Scrapers"compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.327 0.570 0.33 ns Kruskal-Wallis
## 2 Tenericutes Abundance 0.462 0.62 0.46 ns Kruskal-Wallis
## 3 Planctomycetes Abundance 0.0500 0.35 0.05 ns Kruskal-Wallis
## 4 Actinobacteria Abundance 0.623 0.62 0.62 ns Kruskal-Wallis
## 5 Firmicutes Abundance 0.623 0.62 0.62 ns Kruskal-Wallis
## 6 Bacteroidetes Abundance 0.221 0.51 0.22 ns Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.142 0.5 0.14 ns Kruskal-WallisMeans<-compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Phylum group1 group2 p
## 3 Planctomycetes Ostana PDR 0.06619258ggplot(dfShredders, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Shredders(>1%)")
ggplot(dfScrapers, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Scrapers (>1%)")
############
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
dfShredders<-subset(df,FFG =="Shredder")
dfScrapers<-subset(df,FFG == "Scraper")
Means<-compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Family group1 group2 p
## 3 Lachnospiraceae Ostana PDR 0.05182993
## 4 Rikenellaceae Ostana PDR 0.05182993
## 6 Rhodobacteraceae Ostana PDR 0.05182993
## 7 Enterobacteriaceae Ostana PDR 0.05182993print("Scrapers")## [1] "Scrapers"compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteriac~ Abundan~ 0.0143 0.086 0.014 * Kruskal-Wa~
## 2 Aeromonadaceae Abundan~ 0.283 0.37 0.283 ns Kruskal-Wa~
## 3 Rhodobacteraceae Abundan~ 0.142 0.34 0.142 ns Kruskal-Wa~
## 4 Comamonadaceae Abundan~ 0.169 0.34 0.169 ns Kruskal-Wa~
## 5 Methylophilaceae Abundan~ 0.371 0.37 0.371 ns Kruskal-Wa~
## 6 Cytophagaceae Abundan~ 0.371 0.37 0.371 ns Kruskal-Wa~
## 7 Desulfovibriona~ Abundan~ NaN NaN NA ? Kruskal-Wa~
## 8 Rikenellaceae Abundan~ NaN NaN NA ? Kruskal-Wa~
## 9 Ruminococcaceae Abundan~ NaN NaN NA ? Kruskal-Wa~
## 10 Lachnospiraceae Abundan~ NaN NaN NA ? Kruskal-Wa~Means<-compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Family group1 group2 p
## 1 Enterobacteriaceae Ostana PDR 0.01996445ggplot(dfShredders, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Shredders(>1%)")
ggplot(dfScrapers, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Scrapers (>1%)")
Family Level
The top 5 Families across all samples are given in the table below
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
dfStats<-subset(df,FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
df2 <- psmelt(FamilyAll)
df2$Abundance<-df2$Abundance*100
FamilyAll<-ddply(df2, c("Family"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
FamilySorted<-FamilyAll[order(-FamilyAll$mean),]
#FamilySorted[1:5,]
FamilyAll2<-ddply(df2, c("Family","FFG"), summarise, #For Rel abu tabs
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
FamilySorted2<-FamilyAll2[order(-FamilyAll$mean),]
FamilySorted2## Family FFG N mean sd se
## 62 Aeromonadaceae Predator 6 0.000000000 0.00000000 0.000000000
## 56 Acidobacteriaceae Shredder 7 0.004761905 0.01259882 0.004761905
## 126 Brevibacteriaceae Predator 6 0.038888889 0.09525793 0.038888889
## 43 A1-B1 Scraper 9 0.000000000 0.00000000 0.000000000
## 49 Acetobacteraceae Filterer 4 0.050000000 0.05773503 0.028867513
## 133 Burkholderiaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 84 Anaeroplasmataceae Shredder 7 0.090476190 0.08758971 0.033105799
## 16 [Exiguobacteraceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 132 Brucellaceae Shredder 7 0.014285714 0.02622653 0.009912695
## 90 Armatimonadaceae Predator 6 0.155555556 0.15587269 0.063634760
## 135 Burkholderiaceae Scraper 9 0.029629630 0.08888889 0.029629630
## 119 Bifidobacteriaceae Scraper 9 0.003703704 0.01111111 0.003703704
## 149 Chamaesiphonaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 95 auto67_4W Scraper 9 0.000000000 0.00000000 0.000000000
## 33 0319-6G20 Filterer 4 0.000000000 0.00000000 0.000000000
## 39 211ds20 Scraper 9 0.000000000 0.00000000 0.000000000
## 69 Alcaligenaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 147 Caulobacteraceae Scraper 9 0.125925926 0.37777778 0.125925926
## 89 Armatimonadaceae Filterer 4 0.016666667 0.01924501 0.009622504
## 117 Bifidobacteriaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 50 Acetobacteraceae Predator 6 0.016666667 0.02788867 0.011385501
## 114 Beijerinckiaceae Predator 6 0.000000000 0.00000000 0.000000000
## 150 Chamaesiphonaceae Predator 6 0.027777778 0.05340273 0.021801574
## 141 Caldicoprobacteraceae Filterer 4 0.000000000 0.00000000 0.000000000
## 98 Bacillaceae Predator 6 0.011111111 0.02721655 0.011111111
## 78 Anaerobrancaceae Predator 6 0.000000000 0.00000000 0.000000000
## 115 Beijerinckiaceae Scraper 9 0.400000000 0.80484643 0.268282144
## 146 Caulobacteraceae Predator 6 0.277777778 0.13109228 0.053518198
## 113 Beijerinckiaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 128 Brevibacteriaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 54 Acidobacteriaceae Predator 6 0.000000000 0.00000000 0.000000000
## 123 Bradyrhizobiaceae Scraper 9 0.122222222 0.25927249 0.086424162
## 73 Alteromonadaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 129 Brucellaceae Filterer 4 0.008333333 0.01666667 0.008333333
## 86 Anaplasmataceae Predator 6 0.027777778 0.06804138 0.027777778
## 3 [Bryobacteraceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 104 Bacteriovoracaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 64 Aeromonadaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 29 [Weeksellaceae] Filterer 4 0.008333333 0.01666667 0.008333333
## 48 A4b Shredder 7 0.000000000 0.00000000 0.000000000
## 37 211ds20 Filterer 4 0.000000000 0.00000000 0.000000000
## 41 A1-B1 Filterer 4 0.000000000 0.00000000 0.000000000
## 28 [Tissierellaceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 42 A1-B1 Predator 6 0.038888889 0.06116160 0.024969117
## 106 Bacteroidaceae Predator 6 0.000000000 0.00000000 0.000000000
## 125 Brevibacteriaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 61 Aeromonadaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 75 Alteromonadaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 6 [Cerasicoccaceae] Predator 6 0.005555556 0.01360828 0.005555556
## 8 [Cerasicoccaceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 140 C111 Shredder 7 0.019047619 0.01781742 0.006734350
## 25 [Tissierellaceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 107 Bacteroidaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 143 Caldicoprobacteraceae Scraper 9 0.000000000 0.00000000 0.000000000
## 118 Bifidobacteriaceae Predator 6 0.000000000 0.00000000 0.000000000
## 57 Actinomycetaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 124 Bradyrhizobiaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 138 C111 Predator 6 0.150000000 0.16295875 0.066527633
## 68 AK1AB1_02E Shredder 7 0.000000000 0.00000000 0.000000000
## 92 Armatimonadaceae Shredder 7 0.095238095 0.11126973 0.042056004
## 26 [Tissierellaceae] Predator 6 0.000000000 0.00000000 0.000000000
## 35 0319-6G20 Scraper 9 0.000000000 0.00000000 0.000000000
## 23 [Mogibacteriaceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 102 Bacteriovoracaceae Predator 6 0.283333333 0.14414499 0.058846945
## 32 [Weeksellaceae] Shredder 7 0.085714286 0.09786072 0.036987874
## 45 A4b Filterer 4 0.000000000 0.00000000 0.000000000
## 139 C111 Scraper 9 0.040740741 0.12222222 0.040740741
## 13 [Exiguobacteraceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 20 [Fimbriimonadaceae] Shredder 7 0.009523810 0.01626500 0.006147593
## 31 [Weeksellaceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 24 [Mogibacteriaceae] Shredder 7 0.352380952 0.36201961 0.136830553
## 44 A1-B1 Shredder 7 0.042857143 0.07382232 0.027902216
## 58 Actinomycetaceae Predator 6 0.000000000 0.00000000 0.000000000
## 101 Bacteriovoracaceae Filterer 4 0.016666667 0.03333333 0.016666667
## 105 Bacteroidaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 144 Caldicoprobacteraceae Shredder 7 0.028571429 0.04049953 0.015307382
## 19 [Fimbriimonadaceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 91 Armatimonadaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 110 Bdellovibrionaceae Predator 6 0.244444444 0.22377237 0.091354688
## 122 Bradyrhizobiaceae Predator 6 0.005555556 0.01360828 0.005555556
## 53 Acidobacteriaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 38 211ds20 Predator 6 0.005555556 0.01360828 0.005555556
## 2 [Bryobacteraceae] Predator 6 0.083333333 0.11105554 0.045338235
## 5 [Cerasicoccaceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 134 Burkholderiaceae Predator 6 0.000000000 0.00000000 0.000000000
## 77 Anaerobrancaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 21 [Mogibacteriaceae] Filterer 4 0.025000000 0.03191424 0.015957118
## 52 Acetobacteraceae Shredder 7 0.023809524 0.04178554 0.015793451
## 103 Bacteriovoracaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 99 Bacillaceae Scraper 9 0.237037037 0.71111111 0.237037037
## 81 Anaeroplasmataceae Filterer 4 0.000000000 0.00000000 0.000000000
## 142 Caldicoprobacteraceae Predator 6 0.000000000 0.00000000 0.000000000
## 11 [Chthoniobacteraceae] Scraper 9 0.166666667 0.28037673 0.093458910
## 130 Brucellaceae Predator 6 0.066666667 0.06666667 0.027216553
## 1 [Bryobacteraceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 72 Alcaligenaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 120 Bifidobacteriaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 46 A4b Predator 6 0.022222222 0.04036867 0.016480441
## 63 Aeromonadaceae Scraper 9 8.733333333 15.74738955 5.249129851
## 74 Alteromonadaceae Predator 6 0.111111111 0.22673202 0.092562958
## 97 Bacillaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 14 [Exiguobacteraceae] Predator 6 0.005555556 0.01360828 0.005555556
## 85 Anaplasmataceae Filterer 4 0.000000000 0.00000000 0.000000000
## 112 Bdellovibrionaceae Shredder 7 0.028571429 0.05245305 0.019825390
## 79 Anaerobrancaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 100 Bacillaceae Shredder 7 0.009523810 0.02519763 0.009523810
## 18 [Fimbriimonadaceae] Predator 6 0.100000000 0.11737878 0.047919686
## 34 0319-6G20 Predator 6 0.016666667 0.04082483 0.016666667
## 65 AK1AB1_02E Filterer 4 0.008333333 0.01666667 0.008333333
## 36 0319-6G20 Shredder 7 0.000000000 0.00000000 0.000000000
## 121 Bradyrhizobiaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 22 [Mogibacteriaceae] Predator 6 0.000000000 0.00000000 0.000000000
## 83 Anaeroplasmataceae Scraper 9 0.000000000 0.00000000 0.000000000
## 82 Anaeroplasmataceae Predator 6 0.000000000 0.00000000 0.000000000
## 94 auto67_4W Predator 6 0.150000000 0.16158933 0.065968567
## 131 Brucellaceae Scraper 9 4.648148148 6.95218014 2.317393379
## 12 [Chthoniobacteraceae] Shredder 7 0.038095238 0.06214848 0.023489918
## 59 Actinomycetaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 109 Bdellovibrionaceae Filterer 4 0.158333333 0.14240006 0.071200031
## 9 [Chthoniobacteraceae] Filterer 4 0.275000000 0.18534253 0.092671263
## 66 AK1AB1_02E Predator 6 0.000000000 0.00000000 0.000000000
## 67 AK1AB1_02E Scraper 9 0.000000000 0.00000000 0.000000000
## 127 Brevibacteriaceae Scraper 9 0.129629630 0.38888889 0.129629630
## 148 Caulobacteraceae Shredder 7 0.085714286 0.22677868 0.085714286
## 88 Anaplasmataceae Shredder 7 0.000000000 0.00000000 0.000000000
## 145 Caulobacteraceae Filterer 4 0.016666667 0.03333333 0.016666667
## 15 [Exiguobacteraceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 40 211ds20 Shredder 7 0.000000000 0.00000000 0.000000000
## 55 Acidobacteriaceae Scraper 9 0.037037037 0.11111111 0.037037037
## 60 Actinomycetaceae Shredder 7 0.009523810 0.02519763 0.009523810
## 87 Anaplasmataceae Scraper 9 0.000000000 0.00000000 0.000000000
## 93 auto67_4W Filterer 4 0.050000000 0.05773503 0.028867513
## 96 auto67_4W Shredder 7 0.000000000 0.00000000 0.000000000
## 108 Bacteroidaceae Shredder 7 0.004761905 0.01259882 0.004761905
## 111 Bdellovibrionaceae Scraper 9 0.114814815 0.23339946 0.077799821
## 136 Burkholderiaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 137 C111 Filterer 4 0.091666667 0.18333333 0.091666667
## 4 [Bryobacteraceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 7 [Cerasicoccaceae] Scraper 9 0.081481481 0.18189673 0.060632243
## 10 [Chthoniobacteraceae] Predator 6 0.333333333 0.36757463 0.150061716
## 17 [Fimbriimonadaceae] Filterer 4 0.016666667 0.03333333 0.016666667
## 27 [Tissierellaceae] Scraper 9 0.003703704 0.01111111 0.003703704
## 30 [Weeksellaceae] Predator 6 0.300000000 0.35402448 0.144529889
## 47 A4b Scraper 9 0.000000000 0.00000000 0.000000000
## 51 Acetobacteraceae Scraper 9 0.088888889 0.16914819 0.056382731
## 70 Alcaligenaceae Predator 6 0.000000000 0.00000000 0.000000000
## 71 Alcaligenaceae Scraper 9 0.029629630 0.07718024 0.025726748
## 76 Alteromonadaceae Shredder 7 0.042857143 0.11338934 0.042857143
## 80 Anaerobrancaceae Shredder 7 0.161904762 0.30088054 0.113722155
## 116 Beijerinckiaceae Shredder 7 0.000000000 0.00000000 0.000000000cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 2%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)
cdataplot#+ scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(
Trtdata2 <- ddply(df, c("Family", "FFG","Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station)
cdataplot2
########
dfStats<-subset(df, FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
compare_means(Abundance ~ FFG, data = dfStats, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")# Differences between shredder and scrapers only## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteriac~ Abundan~ 0.340 3.80e-1 0.34041 ns Kruskal-Wa~
## 2 Aeromonadaceae Abundan~ 0.0516 6.40e-2 0.05155 ns Kruskal-Wa~
## 3 Rhodobacteraceae Abundan~ 0.958 9.60e-1 0.95779 ns Kruskal-Wa~
## 4 Ruminococcaceae Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 5 Desulfovibriona~ Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 6 Lachnospiraceae Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 7 Rikenellaceae Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 8 Methylophilaceae Abundan~ 0.000369 6.10e-4 0.00037 *** Kruskal-Wa~
## 9 Comamonadaceae Abundan~ 0.000818 1.20e-3 0.00082 *** Kruskal-Wa~
## 10 Cytophagaceae Abundan~ 0.000369 6.10e-4 0.00037 *** Kruskal-Wa~#df<-subset(df, Family!="Aeromonadaceae"&Family!="Rhodobacteraceae"&Family!="Enterobacteriaceae")# Remove groups that were not sig by Kruksal-Wallis
#df <- df[which(df$Family!="Aeromonadaceae")] #& df$Family!="Rhodobacteraceae"& df$Family!="Enterobacteriaceae"
#df
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Family", p.adjust.method = "fdr")
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Family))
for (i in levels(Means$Family)){
Tax<-i
TaxAbundance<-subset(Means,Family==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+5,label=vec))+scale_fill_manual(values=cbPalette)
cdataplot
(Trtdata)## Family FFG N mean sd se
## 1 Aeromonadaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 2 Aeromonadaceae Predator 6 0.000000000 0.00000000 0.000000000
## 3 Aeromonadaceae Scraper 9 8.733333333 15.74738955 5.249129851
## 4 Aeromonadaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 5 Comamonadaceae Filterer 4 12.333333333 8.80151502 4.400757511
## 6 Comamonadaceae Predator 6 6.294444444 2.43715833 0.994965723
## 7 Comamonadaceae Scraper 9 0.362962963 0.46050831 0.153502769
## 8 Comamonadaceae Shredder 7 2.314285714 1.02068553 0.385782867
## 9 Cytophagaceae Filterer 4 1.500000000 1.09273696 0.546368482
## 10 Cytophagaceae Predator 6 14.833333333 8.29572849 3.386716973
## 11 Cytophagaceae Scraper 9 0.011111111 0.03333333 0.011111111
## 12 Cytophagaceae Shredder 7 1.171428571 0.61626971 0.232928057
## 13 Desulfovibrionaceae Filterer 4 10.708333333 5.08821259 2.544106297
## 14 Desulfovibrionaceae Predator 6 0.255555556 0.62598071 0.255555556
## 15 Desulfovibrionaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 16 Desulfovibrionaceae Shredder 7 9.923809524 3.25671468 1.230922446
## 17 Enterobacteriaceae Filterer 4 0.025000000 0.05000000 0.025000000
## 18 Enterobacteriaceae Predator 6 8.122222222 15.38135475 6.279411783
## 19 Enterobacteriaceae Scraper 9 24.911111111 39.25302182 13.084340608
## 20 Enterobacteriaceae Shredder 7 1.985714286 2.21926390 0.838802912
## 21 Lachnospiraceae Filterer 4 6.000000000 5.05253878 2.526269391
## 22 Lachnospiraceae Predator 6 0.011111111 0.02721655 0.011111111
## 23 Lachnospiraceae Scraper 9 0.000000000 0.00000000 0.000000000
## 24 Lachnospiraceae Shredder 7 7.909523810 3.33042730 1.258783201
## 25 Methylophilaceae Filterer 4 2.508333333 1.45178689 0.725893447
## 26 Methylophilaceae Predator 6 3.155555556 3.02571693 1.235243766
## 27 Methylophilaceae Scraper 9 0.011111111 0.03333333 0.011111111
## 28 Methylophilaceae Shredder 7 5.090476190 2.19145768 0.828293148
## 29 Rhodobacteraceae Filterer 4 3.750000000 3.46297881 1.731489404
## 30 Rhodobacteraceae Predator 6 7.327777778 3.02591890 1.235326218
## 31 Rhodobacteraceae Scraper 9 3.977777778 6.42209727 2.140699090
## 32 Rhodobacteraceae Shredder 7 2.647619048 2.44198289 0.922982775
## 33 Rikenellaceae Filterer 4 4.725000000 4.26287809 2.131439046
## 34 Rikenellaceae Predator 6 0.000000000 0.00000000 0.000000000
## 35 Rikenellaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 36 Rikenellaceae Shredder 7 7.076190476 3.38677790 1.280081725
## 37 Ruminococcaceae Filterer 4 5.250000000 3.50518135 1.752590675
## 38 Ruminococcaceae Predator 6 0.005555556 0.01360828 0.005555556
## 39 Ruminococcaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 40 Ruminococcaceae Shredder 7 9.690476190 4.86722654 1.839638714TrtdataPlot<-subset(Trtdata, Family!="Aeromonadaceae"&Family!="Rhodobacteraceae"&Family!="Enterobacteriaceae"&Family!="Cytophagaceae")
(TrtdataPlot)## Family FFG N mean sd se
## 5 Comamonadaceae Filterer 4 12.333333333 8.80151502 4.400757511
## 6 Comamonadaceae Predator 6 6.294444444 2.43715833 0.994965723
## 7 Comamonadaceae Scraper 9 0.362962963 0.46050831 0.153502769
## 8 Comamonadaceae Shredder 7 2.314285714 1.02068553 0.385782867
## 13 Desulfovibrionaceae Filterer 4 10.708333333 5.08821259 2.544106297
## 14 Desulfovibrionaceae Predator 6 0.255555556 0.62598071 0.255555556
## 15 Desulfovibrionaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 16 Desulfovibrionaceae Shredder 7 9.923809524 3.25671468 1.230922446
## 21 Lachnospiraceae Filterer 4 6.000000000 5.05253878 2.526269391
## 22 Lachnospiraceae Predator 6 0.011111111 0.02721655 0.011111111
## 23 Lachnospiraceae Scraper 9 0.000000000 0.00000000 0.000000000
## 24 Lachnospiraceae Shredder 7 7.909523810 3.33042730 1.258783201
## 25 Methylophilaceae Filterer 4 2.508333333 1.45178689 0.725893447
## 26 Methylophilaceae Predator 6 3.155555556 3.02571693 1.235243766
## 27 Methylophilaceae Scraper 9 0.011111111 0.03333333 0.011111111
## 28 Methylophilaceae Shredder 7 5.090476190 2.19145768 0.828293148
## 33 Rikenellaceae Filterer 4 4.725000000 4.26287809 2.131439046
## 34 Rikenellaceae Predator 6 0.000000000 0.00000000 0.000000000
## 35 Rikenellaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 36 Rikenellaceae Shredder 7 7.076190476 3.38677790 1.280081725
## 37 Ruminococcaceae Filterer 4 5.250000000 3.50518135 1.752590675
## 38 Ruminococcaceae Predator 6 0.005555556 0.01360828 0.005555556
## 39 Ruminococcaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 40 Ruminococcaceae Shredder 7 9.690476190 4.86722654 1.839638714#vec<-c("a","b","a","b","a","b","a","b","a","b","a","b")
cdataplot=ggplot(TrtdataPlot, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)#theme(axis.text.x = element_text(angle = 45, hjust = 1))+
cdataplot
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Family", p.adjust.method = "fdr")
#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Family))
#SigList<-length(unique(Trtdata$Family))
#SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
#Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Family", p.adjust.method = "fdr")
#for (i in levels(Means$Family)){
# Tax<-i
# TaxAbundance<-subset(Means,Family==i )
# Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
# difference<-TaxAbundance$p.adj
# names(difference)<-Hyphenated
# Letters<-multcompLetters(difference)
#print(Letters)
# SigList[i]<-Letters
#}
#vec<-unlist(SigList)
#vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
#+geom_text(aes(x=FFG, y=mean+se+5,label=vec))
Means=compare_means(Abundance ~ FFG, data = dfStats,
group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
#keeps
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Family group1 group2 p.adj
## 2 Aeromonadaceae Scraper Shredder 0.07600
## 4 Ruminococcaceae Scraper Shredder 0.00075
## 5 Desulfovibrionaceae Scraper Shredder 0.00075
## 6 Lachnospiraceae Scraper Shredder 0.00075
## 7 Rikenellaceae Scraper Shredder 0.00075
## 8 Methylophilaceae Scraper Shredder 0.00076
## 9 Comamonadaceae Scraper Shredder 0.00140
## 10 Cytophagaceae Scraper Shredder 0.00076Genus level
The top 5 genera across all samples are given in the table below
df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
df2 <- psmelt(GenusAll)
df2$Abundance<-df2$Abundance*100
GenusAll<-ddply(df2, c("Genus"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
GenusSorted<-GenusAll[order(-GenusAll$mean),]
GenusSorted[1:5,]## Genus N mean sd se
## 131 Rhodobacter 26 4.178205 4.430901 0.8689712
## 101 Methylotenera 26 2.488462 2.729264 0.5352528
## 88 Leadbetterella 26 1.992308 5.600399 1.0983287
## 58 Dysgonomonas 26 1.851282 3.369883 0.6608884
## 40 Clostridium 52 1.410897 3.143557 0.4359329GenusAll2<-ddply(df2, c("Genus","FFG"), summarise, #For Rel abu tabs
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
GenusSorted2<-GenusAll2[order(-GenusAll$mean),]
dfStats<-subset(df,FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")Trtdata2 <- ddply(df, c("Genus", "FFG","Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station)
cdataplot2
#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Genus))
#SigList<-length(unique(Trtdata$Genus))
#SigLetters2<-vector(length=NComparisons)
Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Genus", p.adjust.method = "fdr")
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Genus", p.adjust.method = "fdr")
#vec<-unlist(lst)
#for (i in levels(Means$Genus)){
# Tax<-i
# TaxAbundance<-subset(Means,Genus==i )
#print(TaxAbundance)
# Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
# difference<-TaxAbundance$p.adj
# names(difference)<-Hyphenated
# Letters<-multcompLetters(difference)
#print(Letters)
# SigList[i]<-Letters
#}
#vec<-unlist(SigList)
#vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)#+geom_text(aes(x=FFG, y=mean+se+10,label=vec))
cdataplot#+ scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(geom_text(aes(x=FFG, y=mean+se+5,label=vec))
#+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
compare_means(Abundance~FFG, data=dfStats,group.by="Genus",method = "kruskal.test",p.adjust.method="fdr")## # A tibble: 3 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacter Abundance 0.958 0.96 0.95776 ns Kruskal-Wallis
## 2 Clostridium Abundance 0.000239 0.00055 0.00024 *** Kruskal-Wallis
## 3 Methylotenera Abundance 0.000369 0.00055 0.00037 *** Kruskal-WallisMeans=compare_means(Abundance ~ FFG, data = dfStats,
group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 3 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobac~ Abunda~ Scraper Shred~ 1.00e+0 1 1.00000 ns Wilco~
## 2 Clostrid~ Abunda~ Scraper Shred~ 2.99e-4 0.00068 0.00030 *** Wilco~
## 3 Methylot~ Abunda~ Scraper Shred~ 4.57e-4 0.00068 0.00046 *** Wilco~MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats,
group.by = "Genus", p.adjust.method = "fdr")
MeansStation## # A tibble: 3 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacter Abundan~ Ostana PDR 0.0135 0.041 0.014 * Wilcox~
## 2 Clostridium Abundan~ Ostana PDR 0.325 0.49 0.325 ns Wilcox~
## 3 Methyloten~ Abundan~ Ostana PDR 0.866 0.87 0.866 ns Wilcox~# #head(Means)
# keeps <- c("Genus","group1","group2","p.adj","method","p.signif")
# keeps=Means[keeps]
# #keeps
#
#
# test3 <- list('Genus'= keeps$Genus,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
# test3= as.data.frame(test3)
# #test3
# FilteredResults<-test3[!(test3$p.adj>0.05),]
# FilteredResultsPhylumAbuAll
Below is a list of the relative abundance of all taxa at the Phylum level by FFG
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)compare_means(Abundance~FFG, data=df,group.by="Phylum",method = "kruskal.test",p.adjust.method="fdr")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.0147 0.021 0.01466 * Kruskal-Wall~
## 2 Tenericutes Abundance 0.0643 0.075 0.06432 ns Kruskal-Wall~
## 3 Bacteroidetes Abundance 0.000400 0.0014 0.00040 *** Kruskal-Wall~
## 4 Firmicutes Abundance 0.000209 0.0014 0.00021 *** Kruskal-Wall~
## 5 Planctomycetes Abundance 0.0127 0.021 0.01274 * Kruskal-Wall~
## 6 Verrucomicrobia Abundance 0.00184 0.0043 0.00184 ** Kruskal-Wall~
## 7 Actinobacteria Abundance 0.118 0.12 0.11762 ns Kruskal-Wall~Means=compare_means(Abundance ~ FFG, data = df,
group.by = "Phylum", p.adjust.method = "fdr")Multiple comparisons for all taxa grouped by Phylum
NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Phylum))
SigList<-length(unique(Trtdata$Phylum))
SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
#for (i in levels(Means$Phylum)){
# Tax<-i
# TaxAbundance<-subset(Means,Phylum==i )
# print(TaxAbundance)
# Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
# difference<-TaxAbundance$p.adj
# names(difference)<-Hyphenated
# Letters<-multcompLetters(difference)
# print(Letters)
#
#}FamilyAbuAll
Below is a list of the relative abundance of all taxa at the Family level by FFG
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)compare_means(Abundance~FFG, data=df,group.by="Family",method = "kruskal.test",p.adjust.method="fdr")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteria~ Abundan~ 0.0518 5.20e-2 0.05176 ns Kruskal-Wa~
## 2 Aeromonadaceae Abundan~ 0.0365 4.60e-2 0.03654 * Kruskal-Wa~
## 3 Cytophagaceae Abundan~ 0.0000527 1.50e-4 5.3e-05 **** Kruskal-Wa~
## 4 Comamonadaceae Abundan~ 0.0000611 1.50e-4 6.1e-05 **** Kruskal-Wa~
## 5 Rhodobacterace~ Abundan~ 0.0464 5.20e-2 0.04640 * Kruskal-Wa~
## 6 Ruminococcaceae Abundan~ 0.0000556 1.50e-4 5.6e-05 **** Kruskal-Wa~
## 7 Desulfovibrion~ Abundan~ 0.0000686 1.50e-4 6.9e-05 **** Kruskal-Wa~
## 8 Lachnospiraceae Abundan~ 0.0000763 1.50e-4 7.6e-05 **** Kruskal-Wa~
## 9 Rikenellaceae Abundan~ 0.000129 2.10e-4 0.00013 *** Kruskal-Wa~
## 10 Methylophilace~ Abundan~ 0.000339 4.80e-4 0.00034 *** Kruskal-Wa~Means=compare_means(Abundance ~ FFG, data = df,
group.by = "Family", p.adjust.method = "fdr")Multiple comparisons for all taxa grouped by Family
NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Family))
SigList<-length(unique(Trtdata$Family))
SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
for (i in levels(Means$Family)){
Tax<-i
TaxAbundance<-subset(Means,Family==i )
print(TaxAbundance)
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
print(Letters)
}## # A tibble: 3 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aeromonad~ Abunda~ Filter~ Scraper 0.158 0.22 0.15751 ns Wilcox~
## 2 Aeromonad~ Abunda~ Predat~ Scraper 0.0820 0.13 0.08197 ns Wilcox~
## 3 Aeromonad~ Abunda~ Scraper Shredd~ 0.0605 0.097 0.06048 ns Wilcox~
## $Letters
## Filterer Predator Scraper Shredder
## "a" "a" "a" "a"
##
## $LetterMatrix
## a
## Filterer TRUE
## Predator TRUE
## Scraper TRUE
## Shredder TRUE
##
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Comamona~ Abunda~ Filter~ Predat~ 2.40e-1 0.31 0.23953 ns Wilco~
## 2 Comamona~ Abunda~ Filter~ Scraper 6.55e-3 0.016 0.00655 ** Wilco~
## 3 Comamona~ Abunda~ Filter~ Shredd~ 1.07e-2 0.023 0.01073 * Wilco~
## 4 Comamona~ Abunda~ Predat~ Scraper 1.69e-3 0.0073 0.00169 ** Wilco~
## 5 Comamona~ Abunda~ Predat~ Shredd~ 5.28e-3 0.014 0.00528 ** Wilco~
## 6 Comamona~ Abunda~ Scraper Shredd~ 9.89e-4 0.0046 0.00099 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "a" "b" "c"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Cytophag~ Abunda~ Filter~ Predat~ 1.42e-2 0.028 0.01421 * Wilco~
## 2 Cytophag~ Abunda~ Filter~ Scraper 2.08e-3 0.0073 0.00208 ** Wilco~
## 3 Cytophag~ Abunda~ Filter~ Shredd~ 9.25e-1 0.94 0.92472 ns Wilco~
## 4 Cytophag~ Abunda~ Predat~ Scraper 7.06e-4 0.0046 0.00071 *** Wilco~
## 5 Cytophag~ Abunda~ Predat~ Shredd~ 3.41e-3 0.0095 0.00341 ** Wilco~
## 6 Cytophag~ Abunda~ Scraper Shredd~ 4.57e-4 0.0043 0.00046 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "c" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Desulfovi~ Abunda~ Filter~ Preda~ 8.91e-3 0.02 0.00891 ** Wilco~
## 2 Desulfovi~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097 *** Wilco~
## 3 Desulfovi~ Abunda~ Filter~ Shred~ 6.37e-1 0.7 0.63660 ns Wilco~
## 4 Desulfovi~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34 0.27630 ns Wilco~
## 5 Desulfovi~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259 ** Wilco~
## 6 Desulfovi~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobact~ Abunda~ Filter~ Predat~ 0.0304 0.052 0.03038 * Wilco~
## 2 Enterobact~ Abunda~ Filter~ Scraper 0.0189 0.035 0.01892 * Wilco~
## 3 Enterobact~ Abunda~ Filter~ Shredd~ 0.0334 0.055 0.03336 * Wilco~
## 4 Enterobact~ Abunda~ Predat~ Scraper 0.768 0.83 0.76828 ns Wilco~
## 5 Enterobact~ Abunda~ Predat~ Shredd~ 0.886 0.94 0.88625 ns Wilco~
## 6 Enterobact~ Abunda~ Scraper Shredd~ 0.368 0.44 0.36791 ns Wilco~
## Filterer Predator Scraper Shredder
## "a" "ab" "b" "ab"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Lachnospi~ Abunda~ Filter~ Preda~ 1.21e-2 0.025 0.01206 * Wilco~
## 2 Lachnospi~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097 *** Wilco~
## 3 Lachnospi~ Abunda~ Filter~ Shred~ 6.37e-1 0.7 0.63660 ns Wilco~
## 4 Lachnospi~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34 0.27630 ns Wilco~
## 5 Lachnospi~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259 ** Wilco~
## 6 Lachnospi~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Methyloph~ Abunda~ Filter~ Preda~ 9.15e-1 0.94 0.91511 ns Wilco~
## 2 Methyloph~ Abunda~ Filter~ Scrap~ 2.08e-3 0.0073 0.00208 ** Wilco~
## 3 Methyloph~ Abunda~ Filter~ Shred~ 1.56e-1 0.22 0.15638 ns Wilco~
## 4 Methyloph~ Abunda~ Predat~ Scrap~ 8.78e-4 0.0046 0.00088 *** Wilco~
## 5 Methyloph~ Abunda~ Predat~ Shred~ 2.25e-1 0.3 0.22464 ns Wilco~
## 6 Methyloph~ Abunda~ Scraper Shred~ 4.57e-4 0.0043 0.00046 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "a" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacte~ Abunda~ Filter~ Predat~ 0.110 0.17 0.10982 ns Wilco~
## 2 Rhodobacte~ Abunda~ Filter~ Scraper 0.589 0.67 0.58915 ns Wilco~
## 3 Rhodobacte~ Abunda~ Filter~ Shredd~ 0.219 0.3 0.21930 ns Wilco~
## 4 Rhodobacte~ Abunda~ Predat~ Scraper 0.0292 0.051 0.02924 * Wilco~
## 5 Rhodobacte~ Abunda~ Predat~ Shredd~ 0.0184 0.035 0.01842 * Wilco~
## 6 Rhodobacte~ Abunda~ Scraper Shredd~ 1 1 1.00000 ns Wilco~
## Filterer Predator Scraper Shredder
## "ab" "a" "ab" "b"
## # A tibble: 5 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rikenell~ Abunda~ Filter~ Predat~ 2.58e-2 0.047 0.02577 * Wilco~
## 2 Rikenell~ Abunda~ Filter~ Scraper 6.67e-3 0.016 0.00667 ** Wilco~
## 3 Rikenell~ Abunda~ Filter~ Shredd~ 3.95e-1 0.46 0.39509 ns Wilco~
## 4 Rikenell~ Abunda~ Predat~ Shredd~ 2.07e-3 0.0073 0.00207 ** Wilco~
## 5 Rikenell~ Abunda~ Scraper Shredd~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Ruminococ~ Abunda~ Filter~ Preda~ 8.91e-3 0.02 0.00891 ** Wilco~
## 2 Ruminococ~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097 *** Wilco~
## 3 Ruminococ~ Abunda~ Filter~ Shred~ 1.56e-1 0.22 0.15638 ns Wilco~
## 4 Ruminococ~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34 0.27630 ns Wilco~
## 5 Ruminococ~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259 ** Wilco~
## 6 Ruminococ~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"GenusAbuAll
Below is a list of the relative abundance of all taxa at the Genus level by Decomp Stage
GenusSorted2## Genus FFG N mean sd se
## 131 Caldicoprobacter Scraper 9 0.000000000 0.000000000 0.000000000
## 101 Bacteroides Filterer 4 0.000000000 0.000000000 0.000000000
## 88 Asticcacaulis Shredder 7 0.000000000 0.000000000 0.000000000
## 58 Anaerospora Predator 6 0.133333333 0.183787317 0.075030858
## 40 Agrobacterium Shredder 7 0.004761905 0.012598816 0.004761905
## 129 Caldicoprobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 133 Caulobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 69 Aquabacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 117 Bradyrhizobium Filterer 4 0.000000000 0.000000000 0.000000000
## 100 Bacillus Shredder 7 0.009523810 0.025197632 0.009523810
## 59 Anaerospora Scraper 9 0.062962963 0.176733041 0.058911014
## 121 Brevibacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 103 Bacteroides Scraper 9 0.000000000 0.000000000 0.000000000
## 6 02d06 Predator 6 0.000000000 0.000000000 0.000000000
## 70 Aquabacterium Predator 6 0.027777778 0.068041382 0.027777778
## 53 Anaerofustis Filterer 4 0.000000000 0.000000000 0.000000000
## 96 Azospirillum Shredder 7 0.000000000 0.000000000 0.000000000
## 79 Armatimonas Scraper 9 0.000000000 0.000000000 0.000000000
## 142 Chlorobaculum Predator 6 0.011111111 0.027216553 0.011111111
## 13 Achromobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 119 Bradyrhizobium Scraper 9 0.122222222 0.259272486 0.086424162
## 77 Armatimonas Filterer 4 0.008333333 0.016666667 0.008333333
## 128 Burkholderia Shredder 7 0.000000000 0.000000000 0.000000000
## 143 Chlorobaculum Scraper 9 0.000000000 0.000000000 0.000000000
## 111 Bifidobacterium Scraper 9 0.003703704 0.011111111 0.003703704
## 125 Burkholderia Filterer 4 0.000000000 0.000000000 0.000000000
## 164 CM44 Shredder 7 0.071428571 0.073102089 0.027629992
## 145 Christensenella Filterer 4 0.033333333 0.066666667 0.033333333
## 158 Clostridium Predator 12 0.000000000 0.000000000 0.000000000
## 43 Alistipes Scraper 9 0.000000000 0.000000000 0.000000000
## 73 Arenimonas Filterer 4 0.008333333 0.016666667 0.008333333
## 27 Actinomyces Scraper 9 0.000000000 0.000000000 0.000000000
## 82 Arthrobacter Predator 6 0.033333333 0.081649658 0.033333333
## 47 Anaerococcus Scraper 9 0.003703704 0.011111111 0.003703704
## 141 Chlorobaculum Filterer 4 0.000000000 0.000000000 0.000000000
## 15 Achromobacter Scraper 9 0.029629630 0.077180244 0.025726748
## 80 Armatimonas Shredder 7 0.047619048 0.066268653 0.025047197
## 17 Acidovorax Filterer 4 0.000000000 0.000000000 0.000000000
## 11 A17 Scraper 9 0.237037037 0.711111111 0.237037037
## 3 [Clostridium] Scraper 9 0.040740741 0.122222222 0.040740741
## 98 Bacillus Predator 6 0.000000000 0.000000000 0.000000000
## 157 Clostridium Filterer 8 2.845833333 4.455634475 1.575304676
## 159 Clostridium Scraper 18 0.001851852 0.007856742 0.001851852
## 94 Azospirillum Predator 6 0.138888889 0.340206909 0.138888889
## 140 Cellvibrio Shredder 7 0.042857143 0.113389342 0.042857143
## 31 Actinoplanes Scraper 9 0.000000000 0.000000000 0.000000000
## 44 Alistipes Shredder 7 0.319047619 0.264475122 0.099962200
## 81 Arthrobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 20 Acidovorax Shredder 7 0.009523810 0.025197632 0.009523810
## 155 Chthoniobacter Scraper 9 0.096296296 0.233002411 0.077667470
## 54 Anaerofustis Predator 6 0.000000000 0.000000000 0.000000000
## 75 Arenimonas Scraper 9 0.000000000 0.000000000 0.000000000
## 30 Actinoplanes Predator 6 0.005555556 0.013608276 0.005555556
## 60 Anaerospora Shredder 7 0.142857143 0.280683218 0.106088285
## 55 Anaerofustis Scraper 9 0.000000000 0.000000000 0.000000000
## 132 Caldicoprobacter Shredder 7 0.028571429 0.040499526 0.015307382
## 149 Chryseobacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 118 Bradyrhizobium Predator 6 0.005555556 0.013608276 0.005555556
## 99 Bacillus Scraper 9 0.000000000 0.000000000 0.000000000
## 108 Bdellovibrio Shredder 7 0.028571429 0.052453053 0.019825390
## 39 Agrobacterium Scraper 9 0.000000000 0.000000000 0.000000000
## 35 Afifella Scraper 9 0.000000000 0.000000000 0.000000000
## 42 Alistipes Predator 6 0.000000000 0.000000000 0.000000000
## 112 Bifidobacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 85 Asticcacaulis Filterer 4 0.000000000 0.000000000 0.000000000
## 24 Acinetobacter Shredder 7 0.047619048 0.063412649 0.023967728
## 29 Actinoplanes Filterer 4 0.008333333 0.016666667 0.008333333
## 151 Chryseobacterium Scraper 9 0.000000000 0.000000000 0.000000000
## 153 Chthoniobacter Filterer 4 0.025000000 0.031914237 0.015957118
## 137 Cellvibrio Filterer 4 0.000000000 0.000000000 0.000000000
## 90 Azospira Predator 6 0.000000000 0.000000000 0.000000000
## 144 Chlorobaculum Shredder 7 0.000000000 0.000000000 0.000000000
## 56 Anaerofustis Shredder 7 0.014285714 0.037796447 0.014285714
## 126 Burkholderia Predator 6 0.000000000 0.000000000 0.000000000
## 91 Azospira Scraper 9 0.007407407 0.022222222 0.007407407
## 49 Anaerofilum Filterer 4 1.025000000 0.412198262 0.206099131
## 116 Brachybacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 154 Chthoniobacter Predator 6 0.000000000 0.000000000 0.000000000
## 41 Alistipes Filterer 4 0.000000000 0.000000000 0.000000000
## 61 Anaerotruncus Filterer 4 0.000000000 0.000000000 0.000000000
## 114 Brachybacterium Predator 6 0.000000000 0.000000000 0.000000000
## 115 Brachybacterium Scraper 9 0.092592593 0.277777778 0.092592593
## 66 Anaerovorax Predator 6 0.000000000 0.000000000 0.000000000
## 124 Brevibacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 113 Brachybacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 163 CM44 Scraper 9 0.000000000 0.000000000 0.000000000
## 19 Acidovorax Scraper 9 0.000000000 0.000000000 0.000000000
## 92 Azospira Shredder 7 0.000000000 0.000000000 0.000000000
## 161 CM44 Filterer 4 0.000000000 0.000000000 0.000000000
## 21 Acinetobacter Filterer 4 0.383333333 0.589726867 0.294863434
## 38 Agrobacterium Predator 6 0.000000000 0.000000000 0.000000000
## 62 Anaerotruncus Predator 6 0.000000000 0.000000000 0.000000000
## 120 Bradyrhizobium Shredder 7 0.000000000 0.000000000 0.000000000
## 71 Aquabacterium Scraper 9 0.000000000 0.000000000 0.000000000
## 105 Bdellovibrio Filterer 4 0.158333333 0.142400062 0.071200031
## 150 Chryseobacterium Predator 6 0.027777778 0.068041382 0.027777778
## 1 [Clostridium] Filterer 4 0.000000000 0.000000000 0.000000000
## 86 Asticcacaulis Predator 6 0.011111111 0.027216553 0.011111111
## 34 Afifella Predator 6 0.016666667 0.040824829 0.016666667
## 46 Anaerococcus Predator 6 0.000000000 0.000000000 0.000000000
## 109 Bifidobacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 162 CM44 Predator 6 0.005555556 0.013608276 0.005555556
## 16 Achromobacter Shredder 7 0.000000000 0.000000000 0.000000000
## 68 Anaerovorax Shredder 7 0.342857143 0.359379241 0.135832586
## 89 Azospira Filterer 4 0.000000000 0.000000000 0.000000000
## 102 Bacteroides Predator 6 0.000000000 0.000000000 0.000000000
## 123 Brevibacterium Scraper 9 0.129629630 0.388888889 0.129629630
## 4 [Clostridium] Shredder 7 0.000000000 0.000000000 0.000000000
## 32 Actinoplanes Shredder 7 0.000000000 0.000000000 0.000000000
## 48 Anaerococcus Shredder 7 0.000000000 0.000000000 0.000000000
## 33 Afifella Filterer 4 0.000000000 0.000000000 0.000000000
## 134 Caulobacter Predator 6 0.000000000 0.000000000 0.000000000
## 104 Bacteroides Shredder 7 0.004761905 0.012598816 0.004761905
## 18 Acidovorax Predator 6 0.000000000 0.000000000 0.000000000
## 52 Anaerofilum Shredder 7 0.519047619 0.304811504 0.115207919
## 130 Caldicoprobacter Predator 6 0.000000000 0.000000000 0.000000000
## 95 Azospirillum Scraper 9 0.000000000 0.000000000 0.000000000
## 37 Agrobacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 74 Arenimonas Predator 6 0.050000000 0.078173596 0.031914237
## 146 Christensenella Predator 6 0.000000000 0.000000000 0.000000000
## 156 Chthoniobacter Shredder 7 0.004761905 0.012598816 0.004761905
## 9 A17 Filterer 4 0.000000000 0.000000000 0.000000000
## 14 Achromobacter Predator 6 0.000000000 0.000000000 0.000000000
## 67 Anaerovorax Scraper 9 0.000000000 0.000000000 0.000000000
## 87 Asticcacaulis Scraper 9 0.000000000 0.000000000 0.000000000
## 136 Caulobacter Shredder 7 0.047619048 0.125988158 0.047619048
## 160 Clostridium Shredder 14 3.611904762 4.123947307 1.102171279
## 51 Anaerofilum Scraper 9 0.000000000 0.000000000 0.000000000
## 110 Bifidobacterium Predator 6 0.000000000 0.000000000 0.000000000
## 139 Cellvibrio Scraper 9 0.000000000 0.000000000 0.000000000
## 2 [Clostridium] Predator 6 0.000000000 0.000000000 0.000000000
## 5 02d06 Filterer 4 0.000000000 0.000000000 0.000000000
## 7 02d06 Scraper 9 0.000000000 0.000000000 0.000000000
## 8 02d06 Shredder 7 0.009523810 0.025197632 0.009523810
## 22 Acinetobacter Predator 6 1.933333333 1.262097020 0.515248951
## 23 Acinetobacter Scraper 9 0.177777778 0.248886409 0.082962136
## 25 Actinomyces Filterer 4 0.000000000 0.000000000 0.000000000
## 36 Afifella Shredder 7 0.000000000 0.000000000 0.000000000
## 83 Arthrobacter Scraper 9 0.022222222 0.066666667 0.022222222
## 84 Arthrobacter Shredder 7 0.000000000 0.000000000 0.000000000
## 93 Azospirillum Filterer 4 0.000000000 0.000000000 0.000000000
## 106 Bdellovibrio Predator 6 0.244444444 0.223772371 0.091354688
## 107 Bdellovibrio Scraper 9 0.114814815 0.233399462 0.077799821
## 127 Burkholderia Scraper 9 0.029629630 0.088888889 0.029629630
## 147 Christensenella Scraper 9 0.000000000 0.000000000 0.000000000
## 10 A17 Predator 6 0.000000000 0.000000000 0.000000000
## 12 A17 Shredder 7 0.000000000 0.000000000 0.000000000
## 26 Actinomyces Predator 6 0.000000000 0.000000000 0.000000000
## 28 Actinomyces Shredder 7 0.009523810 0.025197632 0.009523810
## 45 Anaerococcus Filterer 4 0.000000000 0.000000000 0.000000000
## 50 Anaerofilum Predator 6 0.000000000 0.000000000 0.000000000
## 57 Anaerospora Filterer 4 0.033333333 0.047140452 0.023570226
## 63 Anaerotruncus Scraper 9 0.000000000 0.000000000 0.000000000
## 64 Anaerotruncus Shredder 7 0.042857143 0.065868235 0.024895853
## 65 Anaerovorax Filterer 4 0.000000000 0.000000000 0.000000000
## 72 Aquabacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 76 Arenimonas Shredder 7 0.014285714 0.037796447 0.014285714
## 78 Armatimonas Predator 6 0.144444444 0.162845075 0.066481224
## 97 Bacillus Filterer 4 0.000000000 0.000000000 0.000000000
## 122 Brevibacterium Predator 6 0.038888889 0.095257934 0.038888889
## 135 Caulobacter Scraper 9 0.000000000 0.000000000 0.000000000
## 138 Cellvibrio Predator 6 0.111111111 0.226732017 0.092562958
## 148 Christensenella Shredder 7 0.000000000 0.000000000 0.000000000
## 152 Chryseobacterium Shredder 7 0.033333333 0.074535599 0.028171808PCoA
Weighted Unifrac
-Ellipses represent 95% CI for the mean of each group
physeqStats<-subset_samples(physeq, FFG!="Predator")
physeqStats<-subset_samples(physeqStats,FFG!="Filterer")
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
Jaccard
-Ellipses represent 95% CI for the mean of each group
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station,scales = "free")#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
Bray-curtis
-Ellipses represent 95% CI for the mean of each group
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
PERMANOVAs
Weighted Unifrac
Homogenieity of Multivariate Dispersions
physeqStats<-subset_samples(physeq, FFG!="Predator")
physeqStats<-subset_samples(physeqStats,FFG!="Filterer")
GPdist=phyloseq::distance(physeqStats, "wunifrac")
beta=betadisper(GPdist, sample_data(physeqStats)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00067157 0.00067157 14.985 999 0.001 ***
## Residuals 14 0.00062744 0.00004482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
FFG
GPdist=phyloseq::distance(physeqStats, "wunifrac")
#MONMDS= ordinate(physeq, "NMDS",GPdist)adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.010177 0.0101767 4.9044 0.23867 0.002 **
## Sampling_station 1 0.003762 0.0037616 1.8128 0.08822 0.082 .
## FFG:Sampling_station 1 0.003800 0.0037998 1.8312 0.08912 0.074 .
## Residuals 12 0.024900 0.0020750 0.58399
## Total 15 0.042638 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.010177 0.0101767 4.389 0.23867 0.001 ***
## Residuals 14 0.032462 0.0023187 0.76133
## Total 15 0.042638 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Jaccard
Homogenieity of Multivariate Dispersions
GPdist=phyloseq::distance(physeq, "jaccard")
beta=betadisper(GPdist, sample_data(physeq)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.048886 0.016295 3.9989 999 0.019 *
## Residuals 22 0.089649 0.004075
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
print("Beta dispersion location")## [1] "Beta dispersion location"beta=betadisper(GPdist, sample_data(physeq)$Sampling_station)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0017764 0.00177642 2.0475 999 0.161
## Residuals 24 0.0208227 0.00086761boxplot(beta)
print("Beta dispersion Shredder v Scraper")## [1] "Beta dispersion Shredder v Scraper"physeqShredderVScraper<-subset_samples(physeq, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqShredderVScraper, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVScraper)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.010007 0.010007 3.9245 999 0.071 .
## Residuals 14 0.035700 0.002550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
print("Beta dispersion Shredder v Predator")## [1] "Beta dispersion Shredder v Predator"physeqShredderVPredator<-subset_samples(physeq, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqShredderVPredator, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVPredator)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0017264 0.0017264 0.9972 999 0.343
## Residuals 11 0.0190437 0.0017312boxplot(beta)
print("Beta dispersion Shredder v Filterer")## [1] "Beta dispersion Shredder v Filterer"physeqShredderVFilterer<-subset_samples(physeq, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqShredderVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVFilterer)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.016325 0.0163249 2.8041 999 0.155
## Residuals 9 0.052396 0.0058218boxplot(beta)
print("Beta dispersion Predator v Scraper")## [1] "Beta dispersion Predator v Scraper"physeqPredatorVScraper<-subset_samples(physeq, FFG=="Predator"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPredatorVScraper, "jaccard")
beta=betadisper(GPdist, sample_data(physeqPredatorVScraper)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.002683 0.0026825 0.936 999 0.347
## Residuals 13 0.037257 0.0028659boxplot(beta)
print("Beta dispersion Predator v Filterer")## [1] "Beta dispersion Predator v Filterer"physeqPredatorVFilterer<-subset_samples(physeq, FFG=="Predator"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqPredatorVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqPredatorVFilterer)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.025561 0.0255605 3.7911 999 0.113
## Residuals 8 0.053937 0.0067422boxplot(beta)
print("Beta dispersion Scraper v Filterer")## [1] "Beta dispersion Scraper v Filterer"physeqScraperVFilterer<-subset_samples(physeq, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqScraperVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqScraperVFilterer)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.047159 0.047159 7.3469 999 0.021 *
## Residuals 11 0.070608 0.006419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
FFG
GPdist=phyloseq::distance(physeqStats, "jaccard")
#MONMDS= ordinate(physeq, "NMDS",GPdist)adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.2666 1.26661 3.7600 0.19102 0.001 ***
## Sampling_station 1 0.6445 0.64455 1.9133 0.09721 0.012 *
## FFG:Sampling_station 1 0.6772 0.67721 2.0103 0.10213 0.013 *
## Residuals 12 4.0424 0.33687 0.60964
## Total 15 6.6308 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.2666 1.26661 3.3058 0.19102 0.001 ***
## Residuals 14 5.3642 0.38315 0.80898
## Total 15 6.6308 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Bray-curtis
Homogenieity of Multivariate Dispersions
GPdist=phyloseq::distance(physeqStats, "bray")
beta=betadisper(GPdist, sample_data(physeqStats)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.025993 0.025993 3.8617 999 0.065 .
## Residuals 14 0.094234 0.006731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
FFG
GPdist=phyloseq::distance(physeqStats, "bray")
#MONMDS= ordinate(physeq, "NMDS",GPdist)adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.6644 1.66440 6.4304 0.27136 0.001 ***
## Sampling_station 1 0.6540 0.65397 2.5266 0.10662 0.007 **
## FFG:Sampling_station 1 0.7092 0.70923 2.7401 0.11563 0.008 **
## Residuals 12 3.1060 0.25884 0.50639
## Total 15 6.1336 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.6644 1.66440 5.2138 0.27136 0.001 ***
## Residuals 14 4.4692 0.31923 0.72864
## Total 15 6.1336 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Jaccard Split by Sampling Station
physeqOstana<-subset_samples(physeq, Sampling_station=="Ostana")
physeqOstanaStats<-subset_samples(physeqOstana,FFG!="Predator")
physeqPDRStats<-subset_samples(physeq, Sampling_station!="Ostana")
#MONMDS= ordinate(physeq, "NMDS",GPdist)Ostana (Predator not included in stats due to N = 2)
GPdist=phyloseq::distance(physeqOstanaStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 1.8065 0.90327 2.6038 0.36653 0.001 ***
## Residuals 9 3.1222 0.34691 0.63347
## Total 11 4.9287 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper Ostana")## [1] "Shredder vs Scraper Ostana"physeqOstanaShredderVScraper<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqOstanaShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.86139 0.86139 2.3315 0.27984 0.027 *
## Residuals 6 2.21670 0.36945 0.72016
## Total 7 3.07809 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Filterer Ostana")## [1] "Shredder vs Filterer Ostana"physeqOstanaShredderVFilterer<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))## Set of permutations < 'minperm'. Generating entire set.##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.85695 0.85695 2.8363 0.36194 0.027 *
## Residuals 5 1.51071 0.30214 0.63806
## Total 6 2.36766 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Filterer Ostana")## [1] "Scraper vs Filterer Ostana"physeqOstanaScraperVFilterer<-subset_samples(physeqOstanaStats, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.9765 0.97652 2.7159 0.27953 0.011 *
## Residuals 7 2.5169 0.35956 0.72047
## Total 8 3.4934 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PDR (Includes Predator, Shredder, Scraper)
GPdist=phyloseq::distance(physeqPDRStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 1.8490 0.92448 2.8964 0.3916 0.001 ***
## Residuals 9 2.8726 0.31918 0.6084
## Total 11 4.7216 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper PDR")## [1] "Shredder vs Scraper PDR"physeqPDRShredderVScraper<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPDRShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.0689 1.06893 3.5129 0.36928 0.029 *
## Residuals 6 1.8257 0.30428 0.63072
## Total 7 2.8946 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Predator PDR")## [1] "Shredder vs Predator PDR"physeqPDRShredderVPredator<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.90936 0.90936 2.874 0.32387 0.022 *
## Residuals 6 1.89844 0.31641 0.67613
## Total 7 2.80779 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Predator PDR")## [1] "Scraper vs Predator PDR"physeqPDRScraperVPredator<-subset_samples(physeqPDRStats, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRScraperVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRScraperVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.79514 0.79514 2.3605 0.28234 0.02 *
## Residuals 6 2.02107 0.33684 0.71766
## Total 7 2.81621 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Random Forest
By Sampling Station
set.seed(200)
GPr = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPr = filter_taxa(GPr, function(x) mean(x) > 1e-5, TRUE)
PhylumAll=tax_glom(GPr, "Phylum")
FamilyAll=tax_glom(GPr,"Family")
GenusAll=tax_glom(GPr,"Genus")
ForestData=GenusAll#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$Sampling_station)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 3.85%
## Confusion matrix:
## Ostana PDR class.error
## Ostana 14 0 0.00000000
## PDR 1 11 0.08333333imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:23, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by Sampling Station")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] de46410fbc643aeaed03d6aa3878d71a X1e994910b44683e96c37da4cd04862a4
## [3] e1a9430f9c3611fc149a9f8e65bf5259 X8766563e48c605e235f4ea2485b02ee3
## [5] X768d500dec54c80d5ef769902dadd792 e6cad09c4fd44cdb1b919a77f3d92429
## [7] X237f7729fcc16a8a864b826f36f307c3 X7b9668ebf6830386ee80c560f7b275a7
## [9] a046edb519c4e3cbdf77ada497c4d743 X5e419280ad42e77b964625f286da0f08
## [11] e32f997a31db6624fc4f67c269121d4c X9eab7a1249300ee6e9757c224201f94d
## [13] b4fda0049757d0552cbf3080550ce23b X2a64f7117c72e249a9c2e5850fd6eb3e
## [15] X147463ced695e1a1dc7c68b5c64c89ed X592e7630b2e58f0bcb2a3f9cccb07ceb
## [17] X887b10e0bdad8789f41fb3c687174051 b672de3b577f55a4350eabdce0f18904
## [19] a1be228b7d7fe14508a59d112f71fa1a X129c4eff55474784a99b0c7309c15c50
## [21] ca23f159cd2907c86ea8d9cf10ae5df9 X0fbb434b42e7241efa6451d8f1b429d0
## [23] ad07e5885874179952e16bf2f7bb57b3
## 165 Levels: de46410fbc643aeaed03d6aa3878d71a ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| b4fda0049757d0552cbf3080550ce23b | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodobacterales | Rhodobacteraceae | Anaerospora | NA | NA | NA | NA | NA | NA | NA | NA |
| e1a9430f9c3611fc149a9f8e65bf5259 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodobacterales | Rhodobacteraceae | Rhodobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| ca23f159cd2907c86ea8d9cf10ae5df9 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Hyphomicrobiaceae | Hyphomicrobium | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| e6cad09c4fd44cdb1b919a77f3d92429 | Bacteria | Bacteroidetes | [Saprospirae] | [Saprospirales] | Chitinophagaceae | Sediminibacterium | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
| a1be228b7d7fe14508a59d112f71fa1a | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Rudanella | NA | NA | NA | NA | NA | NA | NA | NA |
| de46410fbc643aeaed03d6aa3878d71a | Bacteria | Verrucomicrobia | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Luteolibacter | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
Trtdata## Genus Sampling_station N mean sd se
## 1 Anaerospora Ostana 14 0.16666667 0.24564584 0.065651613
## 2 Anaerospora PDR 12 0.01388889 0.03320683 0.009585986
## 3 Emticicia Ostana 14 0.50952381 0.54685247 0.146152469
## 4 Emticicia PDR 12 1.06666667 1.39544714 0.402830892
## 5 Hyphomicrobium Ostana 14 0.33571429 0.45750425 0.122273153
## 6 Hyphomicrobium PDR 12 0.18333333 0.40489430 0.116882916
## 7 Leadbetterella Ostana 14 0.14523810 0.26718128 0.071407201
## 8 Leadbetterella PDR 12 4.14722222 7.86069956 2.269188505
## 9 Luteolibacter Ostana 14 0.41904762 0.48739922 0.130262920
## 10 Luteolibacter PDR 12 0.28333333 0.45979794 0.132732231
## 11 Perlucidibaca Ostana 14 0.34047619 0.45726399 0.122208941
## 12 Perlucidibaca PDR 12 1.93055556 3.20673087 0.925703467
## 13 Rhodobacter Ostana 14 4.97142857 5.19945522 1.389612859
## 14 Rhodobacter PDR 12 3.25277778 3.30687861 0.954613627
## 15 Rhodoferax Ostana 14 1.14761905 1.12376917 0.300339943
## 16 Rhodoferax PDR 12 1.31666667 1.24182173 0.358483055
## 17 Rudanella Ostana 14 0.16428571 0.36408515 0.097305850
## 18 Rudanella PDR 12 0.01944444 0.06735753 0.019444444
## 19 Sediminibacterium Ostana 14 0.05476190 0.08534209 0.022808633
## 20 Sediminibacterium PDR 12 0.26944444 0.45736583 0.132030142cdataplot=ggplot(Trtdata, aes(x=Sampling_station,y=mean))+geom_bar(aes(fill = Sampling_station),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Sampling Station")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
compare_means(Abundance ~ Sampling_station, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterella Abundance 0.250 0.5 0.250 ns Kruskal-Wall~
## 2 Rhodobacter Abundance 0.136 0.34 0.136 ns Kruskal-Wall~
## 3 Perlucidibaca Abundance 0.660 0.73 0.660 ns Kruskal-Wall~
## 4 Rhodoferax Abundance 0.796 0.8 0.796 ns Kruskal-Wall~
## 5 Emticicia Abundance 0.531 0.66 0.531 ns Kruskal-Wall~
## 6 Luteolibacter Abundance 0.0435 0.15 0.044 * Kruskal-Wall~
## 7 Sediminibacterium Abundance 0.304 0.51 0.304 ns Kruskal-Wall~
## 8 Hyphomicrobium Abundance 0.411 0.59 0.411 ns Kruskal-Wall~
## 9 Rudanella Abundance 0.0350 0.15 0.035 * Kruskal-Wall~
## 10 Anaerospora Abundance 0.0276 0.15 0.028 * Kruskal-Wall~ByFFGOstana
GenusAllOstanaStats<-subset_samples(GenusAll, Sampling_station=="Ostana")
GenusAllOstanaStats<-subset_samples(GenusAllOstanaStats, FFG!="Predator")
ForestData=GenusAllOstanaStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 8.33%
## Confusion matrix:
## Filterer Scraper Shredder class.error
## Filterer 3 1 0 0.25
## Scraper 0 5 0 0.00
## Shredder 0 0 3 0.00imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by FFG")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] acccb7cec4d146864bc11d37da55dcd0 a7c725b951d62e6f6814fb8ca64a356e
## [3] ad07e5885874179952e16bf2f7bb57b3 X7403bfa23d688a94b469051d9f13ae5f
## [5] X0fbb434b42e7241efa6451d8f1b429d0 b672de3b577f55a4350eabdce0f18904
## [7] a42594aa39e8c4b1efb497c275e791ba a046edb519c4e3cbdf77ada497c4d743
## [9] X7b9668ebf6830386ee80c560f7b275a7 dc7d0d97c9604c01c99c972b878e09a0
## [11] X8766563e48c605e235f4ea2485b02ee3 X07705f7fffe33ad226eab098040fbb75
## [13] X1036caae06a5e5c605ab8120e5b5b7bc e32f997a31db6624fc4f67c269121d4c
## [15] X43dc593942c5c838cdcb10dac0a39474 d9a11b1b600813533e2f5f0cd7a2eb44
## [17] X4367180b0dc40f77a063355ebab6a4bb ca23f159cd2907c86ea8d9cf10ae5df9
## [19] X855eedb949870fe66f915b9293e0de50
## 165 Levels: acccb7cec4d146864bc11d37da55dcd0 ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| a7c725b951d62e6f6814fb8ca64a356e | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Clostridium | NA | NA | NA | NA | NA | NA | NA | NA |
| d9a11b1b600813533e2f5f0cd7a2eb44 | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Coprococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| dc7d0d97c9604c01c99c972b878e09a0 | Bacteria | Firmicutes | Clostridia | Clostridiales | Ruminococcaceae | Ruminococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| ca23f159cd2907c86ea8d9cf10ae5df9 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Hyphomicrobiaceae | Hyphomicrobium | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| acccb7cec4d146864bc11d37da55dcd0 | Bacteria | Proteobacteria | Betaproteobacteria | Methylophilales | Methylophilaceae | Methylotenera | NA | NA | NA | NA | NA | NA | NA | NA |
| a42594aa39e8c4b1efb497c275e791ba | Bacteria | Deferribacteres | Deferribacteres | Deferribacterales | Deferribacteraceae | Mucispirillum | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterel~ Abundance 0.000133 0.00027 0.00013 *** Kruskal-Wal~
## 2 Clostridium Abundance 0.0000465 0.00023 4.6e-05 **** Kruskal-Wal~
## 3 Perlucidibaca Abundance 0.0000695 0.00023 7.0e-05 **** Kruskal-Wal~
## 4 Methylotenera Abundance 0.000339 0.00048 0.00034 *** Kruskal-Wal~
## 5 Mucispirillum Abundance 0.000268 0.00045 0.00027 *** Kruskal-Wal~
## 6 Rhodoferax Abundance 0.000816 0.001 0.00082 *** Kruskal-Wal~
## 7 Emticicia Abundance 0.000131 0.00027 0.00013 *** Kruskal-Wal~
## 8 Ruminococcus Abundance 0.0000225 0.00023 2.3e-05 **** Kruskal-Wal~
## 9 Coprococcus Abundance 0.00663 0.0074 0.00663 ** Kruskal-Wal~
## 10 Hyphomicrobi~ Abundance 0.140 0.14 0.13988 ns Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 53 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.025 0.01392 * Wilco~
## 2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.015 0.00653 ** Wilco~
## 3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.51 0.44341 ns Wilco~
## 4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0047 0.00043 *** Wilco~
## 5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.013 0.00528 ** Wilco~
## 6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0052 0.00125 ** Wilco~
## 7 Clostrid~ Abunda~ Filter~ Preda~ 5.74e-3 0.014 0.00574 ** Wilco~
## 8 Clostrid~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0047 0.00097 *** Wilco~
## 9 Clostrid~ Abunda~ Filter~ Shred~ 7.77e-1 0.79 0.77681 ns Wilco~
## 10 Clostrid~ Abunda~ Predat~ Shred~ 2.07e-3 0.0061 0.00207 ** Wilco~
## # ... with 43 more rowsMeans=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
Tax<-i
TaxAbundance<-subset(Means,Genus==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%) Forest") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot
ByFFGPDR
GenusAllPDRStats<-subset_samples(GenusAll, Sampling_station!="Ostana")
ForestData=GenusAllPDRStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 0%
## Confusion matrix:
## Predator Scraper Shredder class.error
## Predator 4 0 0 0
## Scraper 0 4 0 0
## Shredder 0 0 4 0imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by FFG")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] a046edb519c4e3cbdf77ada497c4d743 X2b05a6e9e6c580f7fc168a5bbb29de2d
## [3] acccb7cec4d146864bc11d37da55dcd0 ad07e5885874179952e16bf2f7bb57b3
## [5] X237f7729fcc16a8a864b826f36f307c3 X5e419280ad42e77b964625f286da0f08
## [7] X1036caae06a5e5c605ab8120e5b5b7bc X1e994910b44683e96c37da4cd04862a4
## [9] X07705f7fffe33ad226eab098040fbb75 X6b78b6e57cad2b7f9420cee2c46db2aa
## [11] edb114398ced5bf958f4f404493a6642 e32f997a31db6624fc4f67c269121d4c
## [13] a7c725b951d62e6f6814fb8ca64a356e X4367180b0dc40f77a063355ebab6a4bb
## [15] e1a9430f9c3611fc149a9f8e65bf5259 de46410fbc643aeaed03d6aa3878d71a
## [17] b672de3b577f55a4350eabdce0f18904 dc7d0d97c9604c01c99c972b878e09a0
## [19] X0fbb434b42e7241efa6451d8f1b429d0
## 165 Levels: a046edb519c4e3cbdf77ada497c4d743 ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| a7c725b951d62e6f6814fb8ca64a356e | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Clostridium | NA | NA | NA | NA | NA | NA | NA | NA |
| dc7d0d97c9604c01c99c972b878e09a0 | Bacteria | Firmicutes | Clostridia | Clostridiales | Ruminococcaceae | Ruminococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| e1a9430f9c3611fc149a9f8e65bf5259 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodobacterales | Rhodobacteraceae | Rhodobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| acccb7cec4d146864bc11d37da55dcd0 | Bacteria | Proteobacteria | Betaproteobacteria | Methylophilales | Methylophilaceae | Methylotenera | NA | NA | NA | NA | NA | NA | NA | NA |
| edb114398ced5bf958f4f404493a6642 | Bacteria | Bacteroidetes | [Saprospirae] | [Saprospirales] | Saprospiraceae | Haliscomenobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
| de46410fbc643aeaed03d6aa3878d71a | Bacteria | Verrucomicrobia | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Luteolibacter | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
#cdataplot
compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterella Abundan~ 1.33e-4 0.00022 0.00013 *** Kruskal-Wal~
## 2 Rhodobacter Abundan~ 4.63e-2 0.046 0.04634 * Kruskal-Wal~
## 3 Clostridium Abundan~ 4.65e-5 0.00022 4.6e-05 **** Kruskal-Wal~
## 4 Perlucidibaca Abundan~ 6.95e-5 0.00022 7.0e-05 **** Kruskal-Wal~
## 5 Methylotenera Abundan~ 3.39e-4 0.00048 0.00034 *** Kruskal-Wal~
## 6 Rhodoferax Abundan~ 8.16e-4 0.001 0.00082 *** Kruskal-Wal~
## 7 Emticicia Abundan~ 1.31e-4 0.00022 0.00013 *** Kruskal-Wal~
## 8 Ruminococcus Abundan~ 2.25e-5 0.00022 2.3e-05 **** Kruskal-Wal~
## 9 Luteolibacter Abundan~ 1.09e-2 0.012 0.01091 * Kruskal-Wal~
## 10 Haliscomenobac~ Abundan~ 1.16e-4 0.00022 0.00012 *** Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 55 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.025 0.01392 * Wilco~
## 2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.015 0.00653 ** Wilco~
## 3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.51 0.44341 ns Wilco~
## 4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0042 0.00043 *** Wilco~
## 5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.014 0.00528 ** Wilco~
## 6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0057 0.00125 ** Wilco~
## 7 Rhodobac~ Abunda~ Filter~ Preda~ 1.10e-1 0.15 0.10982 ns Wilco~
## 8 Rhodobac~ Abunda~ Filter~ Scrap~ 5.89e-1 0.65 0.58915 ns Wilco~
## 9 Rhodobac~ Abunda~ Filter~ Shred~ 2.18e-1 0.28 0.21825 ns Wilco~
## 10 Rhodobac~ Abunda~ Predat~ Scrap~ 2.92e-2 0.046 0.02924 * Wilco~
## # ... with 45 more rowsMeans=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
Tax<-i
TaxAbundance<-subset(Means,Genus==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%) Alpine Prairie") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot
ByFFGCombined
#GenusAllStats<-subset_samples(GenusAll, FFG!="Predator")
#GenusAllStats<-subset_samples(GenusAllStats, FFG!="Filterer")
GenusAllStats<-GenusAll
ForestData=GenusAllStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 3.85%
## Confusion matrix:
## Filterer Predator Scraper Shredder class.error
## Filterer 3 1 0 0 0.25
## Predator 0 6 0 0 0.00
## Scraper 0 0 9 0 0.00
## Shredder 0 0 0 7 0.00imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by FFG")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] ad07e5885874179952e16bf2f7bb57b3 X0fbb434b42e7241efa6451d8f1b429d0
## [3] b672de3b577f55a4350eabdce0f18904 X43dc593942c5c838cdcb10dac0a39474
## [5] e32f997a31db6624fc4f67c269121d4c acccb7cec4d146864bc11d37da55dcd0
## [7] X4367180b0dc40f77a063355ebab6a4bb a7c725b951d62e6f6814fb8ca64a356e
## [9] a046edb519c4e3cbdf77ada497c4d743 dc7d0d97c9604c01c99c972b878e09a0
## [11] X2b05a6e9e6c580f7fc168a5bbb29de2d a42594aa39e8c4b1efb497c275e791ba
## [13] edb114398ced5bf958f4f404493a6642 X1036caae06a5e5c605ab8120e5b5b7bc
## [15] X7403bfa23d688a94b469051d9f13ae5f X6b78b6e57cad2b7f9420cee2c46db2aa
## [17] ef53cef3b90f7b7b0189b52eaf440bfd X237f7729fcc16a8a864b826f36f307c3
## [19] X07705f7fffe33ad226eab098040fbb75
## 165 Levels: ad07e5885874179952e16bf2f7bb57b3 ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| a7c725b951d62e6f6814fb8ca64a356e | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Clostridium | NA | NA | NA | NA | NA | NA | NA | NA |
| dc7d0d97c9604c01c99c972b878e09a0 | Bacteria | Firmicutes | Clostridia | Clostridiales | Ruminococcaceae | Ruminococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| acccb7cec4d146864bc11d37da55dcd0 | Bacteria | Proteobacteria | Betaproteobacteria | Methylophilales | Methylophilaceae | Methylotenera | NA | NA | NA | NA | NA | NA | NA | NA |
| a42594aa39e8c4b1efb497c275e791ba | Bacteria | Deferribacteres | Deferribacteres | Deferribacterales | Deferribacteraceae | Mucispirillum | NA | NA | NA | NA | NA | NA | NA | NA |
| ef53cef3b90f7b7b0189b52eaf440bfd | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Runella | NA | NA | NA | NA | NA | NA | NA | NA |
| edb114398ced5bf958f4f404493a6642 | Bacteria | Bacteroidetes | [Saprospirae] | [Saprospirales] | Saprospiraceae | Haliscomenobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
Trtdata## Genus FFG N mean sd se
## 1 Clostridium Filterer 4 5.68333333 4.98520032 2.49260016
## 2 Clostridium Predator 6 0.00000000 0.00000000 0.00000000
## 3 Clostridium Scraper 9 0.00000000 0.00000000 0.00000000
## 4 Clostridium Shredder 7 6.82380952 3.53619980 1.33655789
## 5 Emticicia Filterer 4 0.76666667 0.36004115 0.18002057
## 6 Emticicia Predator 6 2.06666667 1.42439071 0.58150507
## 7 Emticicia Scraper 9 0.01111111 0.03333333 0.01111111
## 8 Emticicia Shredder 7 0.62380952 0.47442530 0.17931591
## 9 Haliscomenobacter Filterer 4 0.07500000 0.07391186 0.03695593
## 10 Haliscomenobacter Predator 6 0.13888889 0.12003086 0.04900239
## 11 Haliscomenobacter Scraper 9 0.00000000 0.00000000 0.00000000
## 12 Haliscomenobacter Shredder 7 0.00000000 0.00000000 0.00000000
## 13 Leadbetterella Filterer 4 0.10000000 0.11547005 0.05773503
## 14 Leadbetterella Predator 6 8.38333333 9.64634070 3.93810210
## 15 Leadbetterella Scraper 9 0.00000000 0.00000000 0.00000000
## 16 Leadbetterella Shredder 7 0.15714286 0.14104673 0.05331065
## 17 Methylotenera Filterer 4 2.50833333 1.45178689 0.72589345
## 18 Methylotenera Predator 6 3.15555556 3.02571693 1.23524377
## 19 Methylotenera Scraper 9 0.01111111 0.03333333 0.01111111
## 20 Methylotenera Shredder 7 5.09047619 2.19145768 0.82829315
## 21 Mucispirillum Filterer 4 2.76666667 2.05065482 1.02532741
## 22 Mucispirillum Predator 6 0.01666667 0.04082483 0.01666667
## 23 Mucispirillum Scraper 9 0.00000000 0.00000000 0.00000000
## 24 Mucispirillum Shredder 7 0.81428571 0.57214902 0.21625200
## 25 Perlucidibaca Filterer 4 0.75833333 0.30107831 0.15053915
## 26 Perlucidibaca Predator 6 4.11111111 3.41107130 1.39256403
## 27 Perlucidibaca Scraper 9 0.00000000 0.00000000 0.00000000
## 28 Perlucidibaca Shredder 7 0.03333333 0.05091751 0.01924501
## 29 Rhodoferax Filterer 4 2.19166667 1.59591492 0.79795746
## 30 Rhodoferax Predator 6 2.16111111 0.97260627 0.39706485
## 31 Rhodoferax Scraper 9 0.18518519 0.32236587 0.10745529
## 32 Rhodoferax Shredder 7 1.20952381 0.52200266 0.19729846
## 33 Ruminococcus Filterer 4 0.00000000 0.00000000 0.00000000
## 34 Ruminococcus Predator 6 0.00000000 0.00000000 0.00000000
## 35 Ruminococcus Scraper 9 0.00000000 0.00000000 0.00000000
## 36 Ruminococcus Shredder 7 1.16666667 1.20308246 0.45472243
## 37 Runella Filterer 4 0.27500000 0.44253060 0.22126530
## 38 Runella Predator 6 0.70555556 0.91394546 0.37311667
## 39 Runella Scraper 9 0.00000000 0.00000000 0.00000000
## 40 Runella Shredder 7 0.01428571 0.03779645 0.01428571cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterella Abundan~ 1.33e-4 0.00022 0.00013 *** Kruskal-Wal~
## 2 Clostridium Abundan~ 4.65e-5 0.00022 4.6e-05 **** Kruskal-Wal~
## 3 Perlucidibaca Abundan~ 6.95e-5 0.00022 7.0e-05 **** Kruskal-Wal~
## 4 Methylotenera Abundan~ 3.39e-4 0.00038 0.00034 *** Kruskal-Wal~
## 5 Mucispirillum Abundan~ 2.68e-4 0.00036 0.00027 *** Kruskal-Wal~
## 6 Rhodoferax Abundan~ 8.16e-4 0.00082 0.00082 *** Kruskal-Wal~
## 7 Emticicia Abundan~ 1.31e-4 0.00022 0.00013 *** Kruskal-Wal~
## 8 Ruminococcus Abundan~ 2.25e-5 0.00022 2.3e-05 **** Kruskal-Wal~
## 9 Runella Abundan~ 2.89e-4 0.00036 0.00029 *** Kruskal-Wal~
## 10 Haliscomenobac~ Abundan~ 1.16e-4 0.00022 0.00012 *** Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 55 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.022 0.01392 * Wilco~
## 2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.013 0.00653 ** Wilco~
## 3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.48 0.44341 ns Wilco~
## 4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0036 0.00043 *** Wilco~
## 5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.012 0.00528 ** Wilco~
## 6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0046 0.00125 ** Wilco~
## 7 Clostrid~ Abunda~ Filter~ Preda~ 5.74e-3 0.012 0.00574 ** Wilco~
## 8 Clostrid~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0041 0.00097 *** Wilco~
## 9 Clostrid~ Abunda~ Filter~ Shred~ 7.77e-1 0.79 0.77681 ns Wilco~
## 10 Clostrid~ Abunda~ Predat~ Shred~ 2.07e-3 0.0054 0.00207 ** Wilco~
## # ... with 45 more rowsMeans=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
Tax<-i
TaxAbundance<-subset(Means,Genus==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot
PiCRUST
PiCRUST<-import_biom("C:\\Users\\Joe Receveur\\Downloads\\ItalyInvertPiCRUST2.20.19.biom")
PiCRUST<-merge_phyloseq(PiCRUST,sampdat)
PiCRUST=rarefy_even_depth(PiCRUST, 750000, replace = TRUE, trimOTUs = TRUE, verbose = TRUE,rngseed = TRUE)## `set.seed(TRUE)` was used to initialize repeatable random subsampling.## Please record this for your records so others can reproduce.## Try `set.seed(TRUE); .Random.seed` for the full vector## ...## 487OTUs were removed because they are no longer
## present in any sample after random subsampling## ...PiCRUST by sampling x FFG
GPdist=phyloseq::distance(PiCRUST, "jaccard")
adonis(GPdist ~ FFG*Sampling_station, as(sample_data(PiCRUST), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(PiCRUST), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 3 0.90035 0.300117 5.0869 0.40838 0.001 ***
## Sampling_station 1 0.08310 0.083101 1.4085 0.03769 0.175
## FFG:Sampling_station 2 0.10029 0.050145 0.8499 0.04549 0.641
## Residuals 19 1.12096 0.058998 0.50844
## Total 25 2.20470 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PiCRUST beta div by FFG
print("Shredder vs Scraper")## [1] "Shredder vs Scraper"PiCRUSTShredderVScraper<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(PiCRUSTShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.50616 0.50616 6.4478 0.31533 0.001 ***
## Residuals 14 1.09902 0.07850 0.68467
## Total 15 1.60518 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Filterer")## [1] "Shredder vs Filterer"PiCRUSTShredderVFilterer<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTShredderVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTShredderVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.057191 0.057191 2.5935 0.2237 0.027 *
## Residuals 9 0.198462 0.022051 0.7763
## Total 10 0.255653 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Predator")## [1] "Shredder vs Predator"physeqShredderVPredator<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqShredderVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqShredderVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.33225 0.33225 17.331 0.61173 0.002 **
## Residuals 11 0.21088 0.01917 0.38827
## Total 12 0.54313 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Filterer")## [1] "Scraper vs Filterer"PiCRUSTScraperVFilterer<-subset_samples(PiCRUST, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTScraperVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTScraperVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.29397 0.293967 2.9572 0.21188 0.014 *
## Residuals 11 1.09347 0.099406 0.78812
## Total 12 1.38744 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Predator")## [1] "Scraper vs Predator"physeqScraperVPredator<-subset_samples(PiCRUST, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqScraperVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqScraperVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.3455 0.34550 4.0614 0.23805 0.002 **
## Residuals 13 1.1059 0.08507 0.76195
## Total 14 1.4514 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Predator vs Filterer")## [1] "Predator vs Filterer"PiCRUSTPredatorVFilterer<-subset_samples(PiCRUST, FFG=="Predator"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTPredatorVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTPredatorVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTPredatorVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.13581 0.135811 5.2914 0.39811 0.013 *
## Residuals 8 0.20533 0.025666 0.60189
## Total 9 0.34114 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PiCRUST Jaccard Split by Sampling Station
physeqOstana<-subset_samples(PiCRUST, Sampling_station=="Ostana")
physeqOstanaStats<-subset_samples(physeqOstana,FFG!="Predator")
physeqPDRStats<-subset_samples(PiCRUST, Sampling_station!="Ostana")
#MONMDS= ordinate(physeq, "NMDS",GPdist)Ostana (Predator not included in stats due to N = 2)
GPdist=phyloseq::distance(physeqOstanaStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 0.32232 0.161160 2.2914 0.3374 0.018 *
## Residuals 9 0.63299 0.070333 0.6626
## Total 11 0.95531 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper Ostana")## [1] "Shredder vs Scraper Ostana"physeqOstanaShredderVScraper<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqOstanaShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.20183 0.201826 2.257 0.27334 0.11
## Residuals 6 0.53654 0.089423 0.72666
## Total 7 0.73836 1.00000print("Shredder vs Filterer Ostana")## [1] "Shredder vs Filterer Ostana"physeqOstanaShredderVFilterer<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))## Set of permutations < 'minperm'. Generating entire set.##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.03489 0.034890 1.2951 0.20573 0.294
## Residuals 5 0.13470 0.026939 0.79427
## Total 6 0.16959 1.00000print("Scraper vs Filterer Ostana")## [1] "Scraper vs Filterer Ostana"physeqOstanaScraperVFilterer<-subset_samples(physeqOstanaStats, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.22322 0.223224 2.6273 0.2729 0.024 *
## Residuals 7 0.59475 0.084965 0.7271
## Total 8 0.81798 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PDR (Includes Predator, Shredder, Scraper)
GPdist=phyloseq::distance(physeqPDRStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 0.61094 0.305470 6.0445 0.57324 0.001 ***
## Residuals 9 0.45483 0.050537 0.42676
## Total 11 1.06577 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper PDR")## [1] "Shredder vs Scraper PDR"physeqPDRShredderVScraper<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPDRShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.35923 0.35923 5.2268 0.46556 0.041 *
## Residuals 6 0.41237 0.06873 0.53444
## Total 7 0.77161 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Predator PDR")## [1] "Shredder vs Predator PDR"physeqPDRShredderVPredator<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.27869 0.278694 24.559 0.80366 0.04 *
## Residuals 6 0.06809 0.011348 0.19634
## Total 7 0.34678 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Predator PDR")## [1] "Scraper vs Predator PDR"physeqPDRScraperVPredator<-subset_samples(physeqPDRStats, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRScraperVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRScraperVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.27848 0.278483 3.8931 0.39352 0.028 *
## Residuals 6 0.42920 0.071533 0.60648
## Total 7 0.70768 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PCoA PiCRUST
ord=ordinate(PiCRUST,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(PiCRUST,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse